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Abstract 

The conservation equations for simulating hyper- 
sonic flows in thermal and chemical nonequilibrium 
and details of the associated physical models are pre- 
sented. These details include the curve fits used for 
defining thermodynamic properties of the 11-species 
air model (N, O, N 2 , O 2 , NO, N + , 0 + , N^" , O 2 , 
NO + , and e“), the curve fits for collision cross sec- 
tions, the expressions for transport properties, the 
chemical kinetic models, and the vibrational and 
electronic energy relaxation models. The expres- 
sions are formulated in the context of either a two- 
or three-temperature model. Greater emphasis is 
placed on the two-temperature model, in which it 
is assumed that the translational and rotational en- 
ergy modes are in equilibrium at the translational 
temperature, and the vibrational, electronic, and 
electron translational energy modes are in equilib- 
rium at the vibrational temperature. The eigenvalues 
and eigenvectors associated with the Jacobian of the 
flux vector are also presented in order to accommo- 
date the “upwind” based numerical solutions of the 
complete equation set. 

Introduction 

Future plans for space transportation and explo- 
ration call for mission trajectories with both sus- 
tained and maneuvering hypersonic flight in the 
Earth’s atmosphere at altitudes greater than 70 km 
and velocities greater than 9 km/s (ref. 1). Aero- 
assisted orbital transfer vehicles will use this domain 
in returning from geosynchronous Earth orbit to low 
Earth orbit for rendezvous with Space Station Free- 
dom. Lunar, planetary, and comet sample-return 
missions will utilize the Earth’s upper atmosphere 
for aerobraking as well. Advanced hypersonic air- 
breathing cruise vehicles may ultimately be called 
on to fly through this domain. The trajectories for 
these missions include flow regimes ranging from con- 
tinuum to free molecular. Substantial portions of 
these trajectories, in the transitional regime between 
free molecular and continuum, will carry the vehicle 
through conditions resulting in chemical and ther- 
mal nonequilibrium within the surrounding shock 
layer. Also, chemical nonequilibrium effects can be 
important well into the continuum regime. 

Nonequilibrium processes occur in a flow when 
the time required for a process to accommodate it- 
self to local conditions within some region is of the 
same order as the transit time across the region. If 
the accommodation time is very short compared with 
the transit time, the process is in equilibrium. If 
the accommodation time is very long compared with 
the transit time, the process is frozen. The length 


scale of the region depends on what is being studied. 
Useful length scales include shock standoff distance, 
shock-transition-zone thickness, and boundary-layer 
thickness. The process can be a chemical reaction 
or an exchange of energy among the various modes 
(translational, rotational, vibrational, or electronic) 
of the atoms and molecules in the gas flow. The ac- 
commodation is manifested through collisions among 
the atoms, molecules, ions, and electrons within the 
gas and the accommodation time is determined by 
the frequency with which effective collisions occur. 
The combination of low density in the upper at- 
mosphere (which lowers the collision frequency) and 
high vehicle velocity (which lowers the transit time) 
creates the conditions which make nonequilibrium 
phenomena an important aspect of the shock-layer 
flow. 

Nonequilibrium processes impact the flow envi- 
ronment over a vehicle in several important ways. 
They can significantly alter the shock standoff dis- 
tance and shock shape. The effective isentropic expo- 
nent of the gas is changed, which in turn affects pres- 
sure distributions over expansion and compression 
surfaces of the vehicle. Thermal nonequilibrium, the 
condition wherein the distribution of energy among 
the translational, rotational, vibrational, and elec- 
tronic modes of the gas cannot be described by a 
single temperature, influences the rates at which cer- 
tain chemical reactions can proceed. Translational 
temperature behind the shock is increased, but vi- 
brational and electronic temperatures are decreased 
because of the finite relaxation times for energy 
transfer caused by chemical and thermal relaxation 
processes. The onset of ionization is enhanced, rel- 
ative to the thermal equilibrium state, because of 
the functional dependence of ionizing reactions on 
the translational temperature, whereas dissociation 
is diminished because of the functional dependence of 
dissociative reactions on the vibrational temperature 
(ref. 2). These effects compete in the determina- 
tion of thermal radiation (ref. 3), which may play 
a significant role in the total heating load encoun- 
tered by a large aeroassisted vehicle returning from 
geosynchronous Earth orbit (or beyond) . 

Most of the previous work on computational sim- 
ulation of continuum, nonequilibrium, hypersonic 
flows has concentrated on chemical nonequilibrium. 
These simulations are based on inviscid equations 
(refs. 4 and 5), boundary- layer equations (refs. 6 
to 8), viscous-shock-layer equations (refs. 9 to 11), 
parabolized Navier-Stokes equations (refs. 12 and 
13), and Navier-Stokes equations (refs. 14 to 18). 
All the continuum simulations employ a chemical ki- 
netic model for air which is used to define the pro- 
duction terms in the species continuity equations. 



By contrast, noncontinuum simulations of hypersonic 
flows using a Direct- Simulation Monte Carlo algo- 
rithm intrinsically account for the effects of both 
chemical and thermal nonequilibrium, provided the 
gas components are modeled with the appropriate 
degrees of freedom (refs. 19 to 21). Some of the 
continuum work is focused particularly on determin- 
ing electron number densities in the flow field in or- 
der to predict radio blackout during entry (refs. 22 
to 24). The chemical kinetic models have evolved 
from both theoretical and experimental investiga- 
tions (ref. 25). A review of both chemical and ther- 
mal nonequilibrium effects in nozzles (ref. 26) reflects 
the state of the art in 1967. More recently, quasi one- 
dimensional nozzle flow simulations have been im- 
plemented which include the effects of thermal non- 
equilibrium with a separate vibrational temperature 
for each species (ref. 27). Although no simulations 
are presented, equation sets and relaxation models 
for flows in chemical and thermal nonequilibrium 
were presented by Lee (ref. 2) and, earlier, by Ap- 
pleton and Bray (ref. 28) for a simpler, monatomic, 
ionizing gas. These equation sets include a three- 
temperature model, in which it is assumed that there 
is a single temperature which describes the distribu- 
tion of energy in the translational- rotational modes, 
a second temperature for the vibrational modes, and 
a third temperature for the electronic-free electron 
modes, and a two-temperature model, in which it is 
assumed that the distribution of energy in both the 
vibrational and electronic modes can be described by 
a single temperature. Much of the recent work on hy- 
personic flows over blunt bodies in both thermal and 
chemical nonequilibrium (refs. 29 to 31) employs the 
kinetic and thermal relaxation models summarized 
in the paper by Lee using one- or two-dimensional 
analyses. 

The equations and models in the present pa- 
per are substantially derived from the work of Park 
(ref. 32) and Lee (ref. 2). Similar contributions 
appear earlier in the literature (refs. 33 and 34). 
Some modifications to and citations of other sources 
have been provided; however, this paper is not in- 
tended to justify every element of the model with 
first-principles statistical mechanics arguments. In 
fact, because of the approximate nature of the two- 
temperature model and because of the evolving un- 
derstanding of relaxation processes through theoret- 
ical and experimental research, it is expected that 
some elements of the model will need to be revised. 
Rather, the model is presented with sufficient ex- 
planation to understand the fundamental assump- 
tions and empiricisms involved in its present state 
of evolution. Also, it is presented from the per- 
spective of a computational fluid dynamicist, and 


therefore consideration is given to the way the conser- 
vation equations and physical models should be for- 
mulated within a numerical algorithm. All the phys- 
ical constants and requisite curve fits are included to 
fully define the nonequilibrium, nonradiating, two- 
temperature model. (The effects of radiation are car- 
ried through the equations, but the description of the 
radiation model itself is beyond the scope of this pa- 
per. The effects of a three-temperature model are 
also carried through, but no data are provided for 
electron-vibrational relaxation rates.) 

The conservation equations are formulated within 
the context of the Langley Aerothermodynamic Up- 
wind Relaxation Algorithm (LAURA, refs. 35 and 
36), and sample calculations are made for a represen- 
tative aeroassisted orbital transfer vehicle trajectory 
point. Only a two-temperature model is considered 
in these calculations. Various aspects of the physi- 
cal model are discussed on the basis of these results. 
Comparisons with a Direct-Simulation Monte Carlo 
algorithm (ref. 37), appropriate for rarefied flows, are 
also presented, and some of the problems involved 
with the use of a continuum-based approximation 
scheme in the transitional region between free molec- 
ular and continuum flows are discussed. It should 
be mentioned that the surface slip-boundary condi- 
tions (refs. 18 and 38) should be employed for ana- 
lyzing flows in the transitional region. These condi- 
tions have not been included in the results presented 
herein. 

Symbols 

A Jacobian matrix of f with respect 

to q 

A s k curve fit constant for evaluating 

C v 

A s constant for determining r^ w 

a frozen sound speed, m/s 

a s e curve fit constant for evaluating 

^, e 

a s curve fit constant for evaluating 

a es 

a sr nondimensional parameter for 

collision of species s and r used 
in defining thermal conductivity 

a® curve fit constant for evaluating 

^v,v 

BA curve fit constant for evaluating 

K Cl r 

b s e curve fit constant for evaluating 
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bs 

curve fit constant for evaluating 


&es 

b s v 

curve fit constant for evaluating 

f-iS 

^v,v 

C b ,r 

preexponential term used in 
evaluating backward reaction rate 
coefficient 

Cf >r 

preexponential term used in 
evaluating forward reaction rate 
coefficient 

C P 

specific heat at constant pressure 
for species s, J/kg-K 

Cp,e 

specific heat at constant pres- 
sure for species s for electronic 
enthalpy,. J /kg-K 

C p ,q 

specific heat at constant pressure 
for mixture for energy mode g, 
where q = t,r, v, e, J/kg-K 

rs 

^p,r 

specific heat at constant pressure 
for species s for rotational 
enthalpy, J/kg-K 

c ’p* 

specific heat at constant pressure 
for species s for translational 
enthalpy, J/kg-K 

C P,V 

specific heat at constant pressure 
for species s for vibrational- 
electronic enthalpy, J/kg-K 

^p,v 

specific heat at constant pressure 
for species s for vibrational 
enthalpy, J/kg-K 

c° 

specific heat at constant volume 
for species s, J/kg-K 

C°, e 

specific heat at constant volume 
for species s for electronic energy, 
J/kg-K 

Cv,q 

specific heat at constant volume 
for mixture for energy mode g, 
where q = J/kg-K 

C s 

specific heat at constant volume 
for species s for rotational 
energy, J/kg-K 

CZ,t 

specific heat at constant volume 
for species s for translational 
energy, J/kg-K 

C S v ,tr 

specific heat at constant volume 
for species s for translational- 
rotational energy, J/kg-K 


C v,V 

specific heat at constant volume 
for species s for vibrational- 
electronic energy, J/kg-K 


specific heat at constant volume 
for species s for vibrational 
energy, J/kg-K 


c s 

curve fit constant for evaluating 


mass fraction of species s 

Cs 

curve fit constant for evaluating 


^es 

C s 

average molecular velocity of 
species s, m/s 

C S 

c v 

curve fit constant for evaluating 
n s 

Cl,C 2 

constants for estimating D s 

D s 

effective diffusion coefficient for 
species s, m 2 /s 

D s 

average vibrational energy per 
unit mass of molecule s, which is 


created or destroyed at rate w s , 

J/kg 

D s 

dissociation energy per unit mass 
of molecule s, J/kg or eV 

D a s 

effective ambipolar diffusion 
coefficient for species s, m 2 /s 

Dgr 

binary diffusion coefficient for 
species s and r, m 2 /s 

E 

total energy per unit mass of 
mixture, J/kg 

E b,r 

activation energy for backward 
reaction r, J 

E f,r 

activation energy for forward 
reaction r, J 

e 

mixture energy per unit mass, 

J/kg 

e e 

mixture electronic energy per 
unit mass, J/kg 

e e,s 

electronic energy per unit mass of 
species s, J/kg 

e s 

energy per unit mass of species s, 

J/kg 


energy of formation of species s, 

J/kg 
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&tr 

mass-weighted average of 
translational-rotational energy 
of reactant species (eq. (54a)), 
j/kg 

kf jT 

e v 

mixture vibrational-electronic 
energy per unit mass, J/kg 

L 

e v 

mass-weighted average of 
vibrational-electronic energy of 

1 


reactant species (eq. (54b)), J/kg 

Ixi ly, I>z 

e v 

mixture vibrational energy per 
unit mass, J/kg 

M s 

e v,s 

vibrational energy (enthalpy) per 
unit mass of species s, J/kg 

m 

e* 

vibrational energy (enthalpy) per 
unit mass of species s evaluated 
at temperature T, J/kg 

m s 

** 

vibrational energy (enthalpy) per 
unit mass of species s evaluated 

m x ,m y ,m z 


at temperature T e , J/kg 

N r 

f 

inviscid flux vector relative to 
computational cell wall 

n 

H 

total enthalpy per unit mass of 
mixture, J/kg 

n b,r 

h 

mixture enthalpy per unit mass, 

J/kg 

ft- e,s 

he,s 

electronic enthalpy per unit mass 
of species s, J/kg 


h s 

enthalpy per unit mass of species 
s, J/kg 

71 f,r 

hs,o 

enthalpy of formation of species 
s, J/kg 


h V,s 

vibrational-electronic enthalpy 
per unit mass of species 5, J/kg 

n 3 

i^x 5 n z 

hv,s 

vibrational enthalpy (energy) per 


unit mass of species 5, J/kg 

P 

h 

first ionization energy of species s 
per kg- mole, J/kg-mole 

Pe 

Ps 

Kc,r 

equilibrium constant for 


reaction r 

Q rad 

k 

Boltzmann’s constant, 
1.380622 x 10" 23 J/K 

a 

kb,r 

backward reaction rate coefficient 
for reaction r, units in cgs 

R 


system consistent with number 
of products 

R 


forward reaction rate coefficient 
for reaction r, units in cgs 
system consistent with number 
of reactants 

matrix of column eigenvectors of 

A 

unit vector tangent to computa- 
tional cell wall 

components of 1 in y , and z 

directions, respectively 

molecular weight of species s, 
kg/kg-mole 

unit vector tangent to computa- 
tional cell wall 

mass of species s per particle, kg 

components of m in x , y, and z 
directions, respectively 

number of reactions in chemical 
kinetic model 

unit vector normal to computa- 
tional cell wall 

exponent of temperature in 
preexponential term for backward 
reaction r 

molar rate of production of 
species s per unit volume by 
electron impact ionization, 
kg-mole/m 3 -s 

exponent of temperature in 
preexponential term for forward 
reaction r 

number density of species j, m -3 

components of n in x, y, and z 
directions, respectively 

pressure, Pa 

electron pressure, Pa 

partial pressure due to species s, 
Pa 

radiative energy transfer rate, 
J/m 3 -s 

vector of conversed variables 

matrix of row eigenvectors of A 

universal gas constant, 

8314.3 J /kg- mole- K 
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R-b,r 

backward reaction rate for 
reaction r, kg-mole/m 3 -s 

Rf,r 

forward reaction rate for 
reaction r, kg-mole/m 3 -s 

T 

translational-rotational 
temperature, K 

T d 

rate-controlling temperature for 
dissociation reactions, K 

T e 

electron-electronic excitation 
temperature, K 

Tp 

average temperature defined in 
equation (44c), K 

Tref 

reference temperature for 
thermodynamic relations, K 

■^sh 

post-shock translational- 
rotational temperature, K 

T V 

vibrational-electron-electronic 
excitation temperature, K 

T v 

vibrational temperature, K 

Ty, sh 

post-shock vibrational 
temperature, K 

t 

time, s 

U 

velocity component normal to 
computational cell wall, m/s 

U 

negative vibrational tempera- 
ture of recombined molecules 
(eqs. (44)), K 

u 

velocity component in 
x - direct ion, m/s 

v? 

velocity vector in three- 
dimensional space, 
j = 1 to 3, m/s 

V 

velocity component tangent 
to computational cell wall in 
1- direct ion, m/s 

V 

vibrational coupling factor 
(eq. (44b)) 

Voo 

free- stream velocity, m/s 

V 

velocity component in 
y-direction, m/s 

W 

velocity component tangent 
to computational cell wall in 
m-direction, m/s 


w 

velocity component in 
^-direction, m/s 

w s 

mass rate of production of 
species s per unit volume, 
kg/m 3 -s 

X 3 

positive vector in three- 
dimensional space, 
j — 1 to 3, m 

Vs 

mole fraction of species s 

Z 

nondimensional temperature used 
in evaluation of if c?r 


partition function for energy 
mode q in species s 

z 

normal distance from body, m 


stoichiometric coefficient for 
reactants in reaction r 

P 

dp 

dpE 

Ps,r 

stoichiometric coefficient for 
products in reaction r 

?ir 

J / k § 

Is 

molar concentration of species s, 
kg- mole /m 3 

a® 

modified collision integrals for 
species s and r, m-s 

8 ij 

Kronecker delta 

V 

frozen thermal conductivity for 
translational-rotational energy of 
heavy particles, J/m-s 

Ve 

frozen thermal conductivity 
for electronic energy due to 
collisions between electrons and 
all particles, J/m-s 

V r 

frozen thermal conductivity for 
rotational energy, J/m-s 

m 

frozen thermal conductivity for 
translational energy, J/m-s 

Vv 

frozen thermal conductivity 
for vibrational energy due to 
collisions between molecules and 
all particles, J/m-s 

A 

diagonal matrix of eigenvalues of 

A 

A 4 

mixture viscosity, N-s/m 2 
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Psj 


u es 


P 

Poo 


Ps 

&es 


°s 


< 'Tes > 


< T S > 


t mw 

*s 


reduced molecular weights 
of species s and j, M s Mj/ 

(M s + My) 

effective collision frequency for 
electrons and heavy particles 
in electronic-translational (e-T) 
energy relaxation, 1/s 

mixture density, kg/m 3 

free-stream mixture density, 
kg/m 3 

density of species s, kg/m 3 

effective electron-neutral energy 
exchange collision cross section 
for species s, m 2 

effective cross section for vibra- 
tional relaxation, m 2 

electronic- vibrational (e-V) 
energy relaxation time for 
molecular species s, s 

translational- vibrational (T-V) 
energy relaxation time for 
molecular species s, s 


fy diffusion-corrected number- 

weighted average translational- 
vibrational (T-V) energy relax- 
ation time for mixture, s 


o(l»l) f>( 2 ’ 2 ) 

i4 sr ) iA sr 


dp 

dpey 


collision integrals for species s 
and r, m 2 


Subscripts: 

The general definitions of subscripts are provided 
below. They are provided to augment the complete 
symbol list, which already includes subscript infor- 
mation. Subscripts may have more than one defini- 
tion, in which case the meaning should be clear from 
the context of the expression and by reference to the 
symbol list above. 


b backward rate quantity 

e electronic mode; electrons 

/ forward rate quantity 

p at constant pressure 

q dummy variable for transla- 

tional, rotational, vibrational, or 
electronic states 

r rotational mode; species r; 

reaction r 


translational- vibrational (T-V) 
energy relaxation time for 
species s from correlation of 
Millikan and White (ref. 65), s 


s species s 

sh post-shock condition 

t translational mode 


translational-vibrational (T-V) 
energy relaxation time for 
species s, limiting form at 
high temperatures from Park 
(ref. 29), s 


V vibrational-electronic mode in 

two-temperature model 

v vibrational mode only; 

at constant volume 

x, y, z Cartesian components 


f v number- weighted average oo free stream 

translational-vibrational (T-V) 
energy relaxation time for 
mixture, s 


Conservation Equations 

The modeled system includes 11 species continuity equations, 3 momentum equations, and 
3 energy equations describing the conservation of vibrational, electronic, and total energies. Species 1 
to 5 are the neutral components of air consisting of N, O, N 2 , O 2 , and NO. Species 6 to 10 are the 
ions corresponding to species 1 to 5, in which one electron has been removed. Species 11 are the 
free electrons. The conservation equations for a reacting gas flow in which thermal nonequilibrium 
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is modeled with a three-temperature approximation (i.e., three-energy equations) can be expressed 
as follows. 


Species Conservation: 


d d . d . _ d , . 

rt l ” + aZ l ” u -^ il,D, dZ i ’ )+ ‘ 


1 2 

Mixture Momentum Conservation: 


d i d i j dp d 

dxl dx % 8x3 


( dij dv?\ 2 du k cii 
M ( dxi + dx { ) 3 M dx fc 


Vibrational Energy Conservation: 


d d j d ( dT v \ d ( ^ » n dj/g 

a/'*’ + (x-^J + (/> E 


( e S,s — e v,s) ( e £*s e «,s) , n 

+ E * < - rs > + E/» < Tes > + 

s=mol. 


s=mol. 


s=mol. 

-✓ N 


5 6 

Electron and Electronic Excitation Energy Conservation : 


d d r j ( .1 jdp e d ( dT e \ d ( , n dy s 

di pee + dxl = U J)xl +d^V e d^J + d^[ P 5 ^ ’ S ^ ^ y 


10 10 /_** _ \ 

+ 2*f*(T - T,) £ ^ - £>«,,?, - £ <* r ~ - <3^ 


s=mol. 
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(i) 


( 2 ) 


(3) 


(4) 


Total Energy Conservation: 


d _ c? __ .■ d ( dT dT v dT e \ 


dx3 J 


d 

dx3 


11 

■D 

S — 1 


<9 2/s 


\ d 


^ du l 

du 3 \ 

2 i du k v,- 

J dx3 

u p j 

^dx3 

dx l J 

- -u p — rS* 3 
3 p dx k 


Qrad 


(5) 


Equations (1) to (5) are taken from reference 2 with some minor changes in notation, the 
addition of a source term in the vibrational energy equation to account for vibrational energy lost 
or gained with dissociation or recombination, and simplifications resulting from assumption of a 
zero conduction current (electron velocity equals ion velocity) and zero charge separation (electron 
number density equals ion number density). A review of the assumptions used in deriving these 
equations and a definition of terms for each equation follow. 



Equations (1) to (5) have been written as par- 
tial differential equations. They are generally found 
in the literature in this form. The discretization 
of a partial differential equation can be formulated 
through the use of either finite differences or finite 
volumes. In general, one associates finite differences 
with the differential form of the conservation laws 
and finite volumes with the integral form of the con- 
servation laws; however, the type of approximation 
scheme need not be tied to any particular representa- 
tion of the conservation equations. In the discussion 
which follows, it is convenient to conceptualize phys- 
ical processes as those which occur in a cell (finite 
volume) or which cross cell walls. The terms which 
describe convective and dissipative processes acting 
across cell walls are expressed in conservative form, in 
which the partial derivative of the quantity is taken 
with respect to x 3 or x % with no leading coefficients. 
All other terms are treated as cell-centered sources 
or sinks of mass and energy. 

Species Conservation 

The four terms in equation (1) represent (1) the 
rate of change of mass of species s per unit volume 
in a cell centered at point x 3 , (2) the flux of mass 
of species s convected across cell walls with mixture 
velocity u 3 , (3) the diffusion of species s across cell 
walls, and (4) the mass production rate of species s 
due to chemical reactions. The mixture density p is 
defined by 

11 

p = Y,p s ( 6 ) 

5=1 

The mole fraction y s is defined by 

y, = -lOi/MlL- m 

E (Pt/M k ) 

k~l 

where the summation limit of 11 is the number of 
species in the mixture and M s is the molecular weight 
of species s. The effective diffusion coefficient D s is 
discussed in the section entitled Transport Properties. 
The production term w s is discussed in the section 
entitled Chemical Kinetic Model . The effects of 
thermal and pressure diffusion have been ignored, 
and the binary diffusion approximation has been 
employed. 

Mixture Momentum Conservation 

The four terms in equation (2) represent (1) the 
rate of change of the ith component of momentum 
per unit volume in a cell centered at point x 3 , (2) the 
flux of the ith component of momentum convected 


across cell walls with mixture velocity u 3 , (3) the 
pressure forces acting on cell walls in the i-direction, 
and (4) the viscous forces acting on cell walls in the 
i-direction. The pressure p is defined by 

11 

P=J2Ps ( 8 ) 

S=1 

and the partial pressure of species s is defined by 

Ps = psR'R y M s (9a) 

where s represents an atomic, molecular, or ionic 
species, and 

p s = p s RT e /M s (9b) 

where s represents the free electron species. In both 
equations, R is the universal gas constant. The 
heavy-particle, translational-rotational temperature 
T and the electron temperature T e are discussed 
in the section entitled Thermodynamic Relations. 
The viscosity /i is discussed in the section entitled 
Transport Properties. Bulk viscosity is assumed to 
be equal to zero in the evaluation of shear stresses. 
An approximation of zero charge separation removes 
an electric field forcing function (proportional to 
dp e /dx l ) from the final expression. 

Vibrational Energy Conservation 

The seven terms in equation (3) represent (1) the 
rate of change of vibrational energy per unit volume 
in a cell centered at point x 3 , (2) the flux of vibra- 
tional energy convected across cell walls with mix- 
ture velocity (3) the conduction of vibrational 
energy across cell walls due to vibrational temper- 
ature gradients, (4) the diffusion of vibrational en- 
ergy across cell walls due to molecular concentration 
gradients, (5) the energy exchange (relaxation) be- 
tween vibrational and translational modes due to col- 
lisions within the cell, (6) the energy exchange (relax- 
ation) between vibrational and electronic modes, and 
(7) the vibrational energy lost or gained due to 
molecular depletion (dissociation) or production (re- 
combination) in the cell. The vibrational energy per 
unit mass e v is defined by 

«. = £ ^ no) 

The vibrational energy per unit mass for species s, 
e V)5 , is defined as a function of T v in the section enti- 
tled Thermodynamic Relations. The vibrational tem- 
perature T v is related to the vibrational energy e v in 
the same section. The vibrational thermal conductiv- 
ity r] v is discussed in the section entitled Transport 
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Properties. The vibrational enthalpy for species s, 
h Vi8 , is identical (ref. 39) to the vibrational energy for 
species s, e VjS , and is used herein to maintain consis- 
tent notation with equation (5). The vibrational en- 
ergies of species s at the translational-rotational tem- 
perature e* s and at the electron temperature e** 5 are 
defined as functions of their temperatures in the sec- 
tion entitled Thermodynamic Relations. The charac- 
teristic relaxation times for translational- vibrational 
(T-V) energy exchange < r s > and for electron- 
vibrational (e-V) energy exchange < r es > are pre- 
sented in the section entitled Relaxation Processes. 
The term D s denotes the vibrational energy level rep- 
resentative of those molecules of species s which are 
preferentially created or destroyed (recombined or 
dissociated) because of their high vibrational quan- 
tum numbers. This quantity is defined in the section 
entitled Chemical Kinetic Model . 

The vibrational energy conservation equation was 
derived under the approximation that the number 
density of all vibrationally excited molecules has 
a Boltzmann distribution characterized by a sin- 
gle temperature T v . The expressions for T-V and 
e-V energy exchange, in which the relaxation rates 
are linearly proportional to the energy difference, 
are based on assumptions of harmonic oscillators 
(ref. 40). Anharmonicity has little effect on the re- 
laxation rate near equilibrium, but the vibrational 
de-excitation rate is enhanced sufficiently far from 
equilibrium (ref. 41). This effect is not included here. 


Electron and Electronic Excitation Energy 
Conservation 


The nine terms in equation (4) represent 
(1) the rate of change of electronic energy per unit 
volume in a cell centered at point x 3 , (2) the flux of 
electronic enthalpy convected across cell walls with 
mixture velocity u 3 , (3) the work done on electrons 
by an electric field induced by the electron pres- 
sure gradient, (4) the conduction of electronic energy 
across cell walls due to the electron temperature gra- 
dient, (5) the diffusion of electronic energy due to 
concentration gradients, (6) the energy exchange due 
to elastic collisions between electrons and heavy par- 
ticles, (7) the energy loss due to electron impact ion- 
ization, (8) the energy exchange (relaxation) due to 
inelastic collisions between electrons and molecules 
in the cell (corresponds to term 6 of eq. (3) for vi- 
brational energy conservation), and (9) rate of energy 
loss due to radiation caused by electronic transitions. 
The electronic energy per unit mass e e contains con- 
tributions from the electronic energy levels of all the 


species and is defined by 
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Ps e e,s 

P 


( 11 ) 


The electronic energy per unit mass of species s, e e<s , 
is defined as a function of T e in the section entitled 
Thermodynamic Relations. The electron tempera- 
ture T e is related to electronic energy e e in that same 
section. The electronic thermal conductivity ry e is 
defined in the section entitled Transport Properties. 
The electronic enthalpy per unit mass of species s, 
h e , s , is identical to e e , s for all species except free 
electrons. In the case of free electrons this quantity 
is defined by 


he,e 


e e,e + 


RT e 

M e , 


( 12 ) 


The effective collision frequency of electrons with 
heavy particles v es is defined in the section entitled 
Relaxation Processes. The molar rate of production 
of species s by electron impact ionization h e ^ s is de- 
fined in the section entitled Chemical Kinetic Model. 
The energy per unit mole lost by a free electron in 
producing species s through electron impact ioniza- 
tion I s is also defined in this section. The radiant en- 
ergy transfer rate due to electronic transitions Q ra d 
is not considered in the present model but is impor- 
tant under some conditions. Consequently, its effects 
should be coupled in subsequent analyses in a man- 
ner similar to that presented in references 42 to 45. 

The electron-electronic excitation energy conser- 
vation equation is derived under the approximation 
that the electronically excited states of all atoms and 
molecules and the translational energies of free elec- 
trons can be characterized by a Maxwell-Boltzmann 
distribution at temperature T e . Furthermore, it is 
approximated that, in the absence of an externally 
applied electric or magnetic field (planetary mag- 
netic field effects are ignored), the charge separation 
in a partially ionized gas is very small because of the 
linking of electron and ion diffusion (ambipolar diffu- 
sion). Also, there is no conduction current in the flow 
field under consideration, and the electron velocity is 
equal to the ion velocity. The inertial terms in the 
electron momentum equation are neglected because 
of the electron’s small mass and viscous stress terms 
due to electrons are assumed to be negligible (refs. 2 
and 28), both of which lead to the approximation 
that the electric field is proportional to the electron 
pressure gradient as used in equation (4). Without 
the assumption of negligible charge separation and 
conduction current, the electron momentum equa- 
tion and the appropriate electrodynamic field terms 
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would have to be used to solve for the electron ve- 
locity vector. The present model cannot account for 
plasma- dynamic effects. These are not expected to 
be important in a continuum, forebody flow field, 
where the estimated magnetic pressures are orders of 
magnitude smaller than the fluid pressures in aero- 
assisted orbital transfer vehicle (AOTV) applica- 
tions. The importance of these effects in the base flow 
region, where electrical conductivity may be higher, 
is not known. 

Total Energy Conservation 

The six terms in equation (5) represent (1) the 
rate of change of total energy per unit volume in a 
cell centered at point x 3 , (2) the flux of total enthalpy 
convected across cell walls with mixture velocity u 3 , 
(3) the conduction of energy across cell walls due to 
temperature gradients, (4) the diffusion of enthalpy 
across cell walls due to concentration gradients, 
(5) the work done by shear forces, and (6) the rate 
of energy loss due to radiation caused by electronic 


transitions. The total energy E is defined by 




irtr 


2 


+ £ 


5=1 


Ps e s 

P 


( 13 ) 


where the energy per unit mass of species s, e 5 , 
is defined in the section entitled Thermodynamic 
Relations. The total enthalpy H is defined by 

H = E+ V - (14) 

P 

The frozen thermal conductivity of heavy particles rj 
is that part of the conductivity arising from collisions 
in which exchanges of translational and rotational 
energy occur. It is defined in the section entitled 
Transport Properties. All other terms in equation (5) 
have been discussed previously. An approximation of 
zero charge separation has been used to simplify this 
equation. 


Two-Temperature Model 

The model developed to this point assumes that the partitioning of energy among the translational, 
rotational, vibrational, and electronic modes in all 11 species can be described adequately by 3 temperatures. 
The relation between these energies and temperatures is discussed in the section entitled Thermodynamic 
Relations. What constitutes an adequate description of the energy distribution is subjective. For the purposes 
of simulating flow fields over hypersonic vehicles, an adequate model is one which allows accurate prediction 
of the aerodynamic coefficients and of both the convective and radiative heating of the vehicle surface. The 
adequacy of models used for the simulation of scramjet engines or gas-dynamic lasers is likely to be judged 
by other requirements. More detailed thermal models can be constructed in which the vibrational energy 
for each molecular species is modeled by its own vibrational temperature (refs. 27 and 31). This treatment 
requires an additional conservation law for the vibrational energy of each molecular species. Candler and 
MacCormack (ref. 31) show relatively small differences between vibrational temperatures of N 2 and O 2 for 
flow over an axisymmetric, blunted cone. Matsuzaki and Hirabayashi (ref. 27) show larger differences among 
vibrational temperatures for expanding flow through a nozzle and for flow behind a normal shock, though 
these temperature distributions are qualitatively similar to each other and are significantly different than the 
translational temperature. These results enhance the credibility of the simpler thermal models, in which the 
vibrational energies of all species are described by a single vibrational temperature. Furthermore, the accuracy 
of a multivibrational temperature model is limited by the accuracy of the available relaxation time data required 
to describe the energy exchange due to collisions among the species. 

Another approach to be considered, particularly in light of uncertainties in some physical parameters and 
the complexities of three-dimensional flow simulation, is to reduce the thermal model to a two-temperature 
system. A justification for a two-temperature model presented by Park (ref. 32) is based on (1) the rapid energy 
transfer between the translational mode of free electrons and the vibrational mode of molecular nitrogen and 
(2) the rapid equilibration of the low-lying electronic states of heavy particles with the ground electronic 
state at the electronic temperature. This model makes the approximation that one temperature, T, describes 
the distribution of heavy-particle translational and rotational energies and that a second temperature, Ty , 
describes the distribution of vibrational, electronic, and electron translational energies. Note that subscript 
V denotes both the vibrational and electronic modes modeled together, whereas subscript v denotes only the 
vibrational modes and subscript e denotes only the electronic modes. While the approximation may be invalid 
in the viscous boundary layer adjacent to the wall (where the physically correct boundary conditions on the 
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vibrational and electronic translational energies are inconsistent with a single temperature (refs. 29 and 46)), 
it allows for a computationally more tractable formulation of reacting flows with thermal nonequilibrium. 

The two-temperature model is obtained by combining equations (3) and (4) for vibrational energy conser- 
vation and electronic energy conservation into a single relation for vibrational-electronic energy conservation, 
where 

T v = T e = T v (15) 

The vibrational- electronic energy conservation equation can now be expressed. 

Vibrational-Electronic Energy Conservation: 


d d , dv? d 

M pev + dZ pevu p 'dJ + d3 


{Vv + Ve) 


dTy 


dx3 


d ( ^ dys ' 

+ ^[ P ^ hv ’ sDs dZ, 

\ s=l / 


+ 


/ sje \ Q 10 10 

V P, { -^^+2pelR(T-T v )J2^-E<A+ E _*A"0»d (16) 
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The vibrational-electronic energy per unit mass, ey, and the vibrational-electronic enthalpy per unit mass for 
species s, hy, s , can be expressed by 

(17) 


ey = e v + e e 
hy,s — hv,s + he,s 


(18) 


Note that terms 2 and 3 of equation (4), dealing with electron pressure, have been combined into a single term 
(term 3) in equation (16) above. The removal of the electron pressure from the convective term simplifies the 
expressions for eigenvalues and eigenvectors of the Jacobian of the flux vector. These quantities are important 
in upwind formulations of the governing equations and are presented in the section entitled Upwind Formulation 
of the Flux Vector. 


Thermodynamic Relations 


where 


If a gas is in thermal equilibrium (i.e., the parti- 
tioning of energy for all modes can be described by a 
Maxwell-Boltzmann distribution at a single temper- 
ature T), then the energy per unit mass of species s 
in the gas can be expressed as 


h s ,o — e s,o + m ^ref 


Cp—Cy + YJ- 


R_ 

M s 


( 21 ) 

( 22 ) 


e s = 



CydT + e S)0 


(19) 


where T re f is a reference temperature, generally taken 
as 298.16 K, C* is the specific heat at constant vol- 
ume for species s, and e Sj0 is the energy of formation 
of species s at temperature T re f. The enthalpy per 
unit mass of species s is similarly expressed by 


hg 



CpdT + h St0 


( 20 ) 


Values for h St0 , M s , and other physical constants 
for the 11 species considered herein are presented 
in table I. Tabulated values and polynomial curve 
fits for Cp and h s , under the assumption of thermal 
equilibrium, are readily available (refs. 47 to 51). The 
curve fits employed for the 11 species considered in 
this report are presented at the end of this section. 

In the general case of thermal nonequilibrium, e s 
and h s are functions of several temperatures, de- 
pending on the number of parameters required to 
adequately describe the partitioning of energy in the 
gas. In the three-temperature model described by 
equations (3) to (5), for example, it is assumed that 
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temperature T describes the translational and rota- 
tional energy modes of heavy particles, T v describes 
the vibrational modes of all molecules, and T e de- 
scribes the electronic excitation and free electron 
translational modes. Consequently, it is necessary 
to establish how e s and h s are functions of these 
temperatures. 

The partition function provides the mechanism 
for establishing these relationships under the as- 
sumptions that there exists a Maxwell- Boltzmann en- 
ergy distribution at the temperature for each mode 
(i.e., translational, rotational, etc.) and that there is 
no coupling of energy levels between modes. Then, 
following the standard methods of statistical mechan- 
ics (refs. 40 and 52), one can write 


_ R d_ ( 2 djn Z^ \ 

V ' q M s dTq \ q dT q ) 


(23) 


where q is a dummy index for the particular energy 
storage mode and is the partition function for 
species s in that mode. For the temperature range of 
interest (200 < T < 50000), both the translational 
and rotational modes are assumed to be fully excited, 
and the specific heat capacity for these modes reduces 
to 


is a coupling between the electronic and vibrational- 
rotational modes in diatomic molecules as well 
because the interatomic forces change when an elec- 
tron leaves the ground electronic state. In this case, 
the partition functions for vibrational and rotational 
energies are taken in the ground electronic state, and 
the electronic energy partition function includes the 
correction to vibration and rotation due to electronic 
transitions evaluated at temperature T e (refs. 52 
to 54). 

Compilations of physical constants for evaluating 
partition functions can be found in references 55 and 
56. Balakrishnan (ref. 53) evaluated the partition 
functions in this manner and generated curve fits for 
the vibrational and electronic heat capacities in the 
following form: 


Cy iV - + K T v + (Vibrational) (26) 

C* e = (^) (aI + 6|T e + i|) (Electronic) (27) 


where the leading factor is a conversion from 
cal/g-mole-K to J/kg-K and the constants a, 6, and 
c are presented in tables 3 and 4 of reference 53. For 
internal degrees of freedom, 


s _ 3 R 
V ' 1 2 M s 

for the translational modes and 


__ 

°v,r — 


R__ 

M's 


(24) 


(25) 


for the rotational modes. 

The partition function for vibrational energy in a 
diatomic molecule is generally derived under the as- 
sumption of a harmonic oscillator, valid for low vibra- 
tional energies (low vibrational quantum numbers), 
with an anharmonic correction required at large vi- 
brational energies. The anharmonic correction is due 
to the effects of interatomic forces on the potential 
energy curve for vibrational energy and due to the 
coupling of rotational and vibrational energies caused 
by a change in the moment of inertia of a diatomic 
molecule with increasing vibrational quantum num- 
bers (ref. 52) . This last effect technically violates the 
assumption of no coupling between modes. The en- 
tire correction is included as part of the vibrational 
energy partition function evaluated at temperature 
T v . 

The partition function for electronic energies is 
obtained by summing over the observed energy level 
data for the atoms and molecules in the gas. There 


C a p , q = C* , (28a) 

where q — r,v, or e and s represents atoms or 
molecules. The electron heat capacities are expressed 
by 


e _ 5 R 
p ’ e " 2 M e 

C e — ^ 
V ’ e 2 M e 


(28b) 

(28c) 


The evaluations of specific heats and enthalpies 
are much simpler in the two-temperature model. The 
curve fits that are available for the enthalpies and 
heat capacities of species as a function of tempera- 
ture are valid only under conditions of thermal equi- 
librium. They assume that a single temperature T 
(where T — T v — T e ) describes the partition of en- 
ergy among all the modes, and it is that tempera- 
ture which must be used in the curve fits. However, 
in the two- temperature model, one can take advan- 
tage of the fact that the translational and rotational 
energy modes are assumed to be fully excited and 
therefore the heat capacities for these modes are inde- 
pendent of temperature. The vibrational-electronic 
heat capacity for species s can be evaluated by uti- 
lizing the curve fit for total heat capacity evaluated 
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at temperature Ty and subtracting out the constant 
contribution from the translational and rotational 
heat capacities. This strategy is employed in equa- 
tions (29) and (30). In like manner, the contribu- 
tion to enthalpy from the translational and rotational 
modes is linear with temperature. Consequently, the 
vibrational-electronic enthalpy for species s can be 
evaluated by utilizing the curve fit for specific en- 
thalpy evaluated at temperature Ty and subtract- 
ing out both the contribution from the translational 
and rotational enthalpy evaluated at Ty and the en- 
thalpy of formation. The correct specific enthalpy 
can then be recovered by adding the contribution 
of translational and rotational enthalpy evaluated 
at temperature T and the enthalpy of formation to 
the vibrational-electronic enthalpy. This strategy is 
used in the derivation of equations (35) and (36). 
The vibrational-electronic heat capacity can now be 
evaluated as follows: 

C s v ,v(T V ) = C S V (T V ) - <7* t - C* r (29) 

where __ 

C?,(T v ) = C;(T v )-~ (30) 

Curve fits and tabulated values for C 3 (T), as 
stated earlier, can be found in references 47 to 51. 
The use of curve fits significantly reduces the com- 
plexity and expense of computing the original func- 
tions. The original expressions for C 3 were obtained 
from either the partition function method outlined 
above or a virial coefficient method. Either source is 
suitable for the two-temperature model. 1 * * * The curve 
fits presented below should not be considered recom- 
mended data. They have been employed because of 
their accessibility and use in other codes or because 
they are the only known fits to data in a given tem- 
perature range. 

The curve fits for C 3 (T) in the present model are 
of the form 

C7 p( T ) = -^E^ Tfc “ 1 ( 31 ) 

3 k = l 

Values for A s k are presented in table II along with 
the original sources. The constants for the two high- 
est temperature ranges in table II are previously 


unpublished data from reference 18. Values of A| 
are linearly averaged across the curve fit boundaries 
(i.e., 950 < T < 1050,5900 < T < 6100,14900 < 
T < 15100, and 24 900 < T < 25100) to ensure 
smooth variation of thermodynamic properties over 
the entire temperature range. 

The evaluation of h s with the three-temperature 
model is obtained by integrating the various heat 
capacities over the appropriate limits and adding the 
heat of formation. Thus, 

hs = [ T (C 8 p>t + Cp r )dT' + F V Cp V dT' 

JT le{ V ’ * JT te f 

+ [ ^ CpgdT 1 + h s , 0 (32) 

Jt tg{ 

where T f is the dummy variable of integration and 
all other terms have been defined. The vibrational 
and electronic energies of species s can be written 
individually as 



C 3 dT' 


e 


e,s 



C 3 dT f 


(33) 

(34) 


The integration of the curve fits for C 3 V and C 3 e 
from equations (26) and (27) is trivial and completes 
the definition. 

In the two-temperature model, obtain hy s and 
h s from equations (29) to (31). Thus, 


h V A T v) = MT«) - (C* t + C*, r )(7V - T re f) - h a ,o (35) 
h s {T,T v ) = h Via (T v ) + (C« t + C* r )(T - T ref ) + h s ,o (36) 

Curve fits for h s (Ty ), which include the heat of 
formation, can be written with equation (31). The 
expression takes the form 


1 The partition function can be difficult to specify for di- 

atomic molecules at vibrational energies near the dissocia- 

tion energy, and to accurately define the anharmonic correc- 

tion near these levels is not a trivial matter. For example, 

Balakrishnan (ref. 53) noted that his partition function eval- 
uation of Cp 2 deviates from the virial coefficient approach 
of Browne (ref. 51) at temperatures greater than 12000 K. 
Research is ongoing to improve these formulations. 


p / ^ ASrpk \ 

(37) 

The constant A| is also provided in table II. The 
value of T re f for these curve fits is 298.16 K. 
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Equations (19) to (37) were developed for a single 
species. The mixture relations are expressed as 


Cp,, = \'£k> c ™ 

(38) 

Cv,q = - PsCv,q 

(39) 

h = — y] ps^s 

(40) 


Equations (19) to (40) and the data in tables I 
and II complete all required thermodynamic relations 
for a two-temperature model. The addition of data 
available in tables 3 and 4 of reference 53 permits 
specification of all thermodynamic relations required 
for a three-temperature model. 

Chemical Kinetic Model 

The mass rate of production of species s per unit 
volume is expressed as 

N r 

W s — M s ^^(/?s,r “ a s,r){Rf,r ~ ^6,r) (41) 

r=l 

where N r is the number of reactions, a Si r and /? Sjr 
are respectively the stoichiometric coefficients for 
reactants and products in the r reaction, and Rj, r 
and R br are respectively the forward and backward 
reaction rates for the r reaction. These rates are 
defined by 


provide these expressions as a function of a single 
temperature. However, under the low-density and 
high-energy flow conditions of interest herein (where 
thermal equilibrium may not be assumed) , the char- 
acteristic chemical time scale for dissociative reac- 
tions is comparable to the characteristic time for 
vibrational relaxation, a condition suggesting a cou- 
pling between the vibrational and chemical processes. 
Models for such chemical-vibrational coupling are 
considered below. 

Chemical-Vibrational Coupling 

Two types of chemical- vibrational coupling have 
been suggested in the literature. Under the first cou- 
pling model, known as preferential dissociation, the 
dissociation of molecular species is obtained more 
easily when the molecules are vibrationally excited. 
Accordingly, the molecules in the higher vibrational 
states are assumed to be preferentially dissociated. 
Molecules in the lower vibrationally excited states 
must “ladder climb” to the higher vibrationally ex- 
cited states before dissociation can occur. However, 
this model may not be valid at very high veloci- 
ties. Under highly energetic conditions the ladder- 
climbing process may not be as significant and a 
second model, based on nonpreferential dissociation, 
may be more realistic. Both models are discussed 
below. 

Preferential Dissociation and Recombination 

Treanor-Marrone model (refs. 57 and 58). In 

this model, the effect of vibrational relaxation on 
dissociation is included through the relation 


R fjV = 1000 


R bi r = 1000 
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k f,r n(°- 0 °WM s )“- 

s = 1 
11 

h,r H(0Mlp s /M s )^ 


8=1 


(42) 

(43) 


where kj, r and fc& r are respectively the forward and 
backward reaction rate coefficients. Reaction rate co- 
efficient data are generally provided in cgs units in 
the literature. The term in brackets is in cgs units. 
The factors 1000 and 0.001 are required in the con- 
version from cgs units to mks units. (Because most 
of the data in the literature for reaction rate coeffi- 
cients are in cgs units, we have retained this prac- 
tice in the formulation of equations (42) and (43).) 
A chemical kinetic model is defined when a set of 
N r reactions is provided with the appropriate expres- 
sions for the forward and backward rate coefficients. 
Most of the sources for reaction rate coefficient data 
have assumed thermal equilibrium and, consequently, 


kf, r = k* f>r V(T,T v ,U) (44a) 


where the vibrational coupling factor V is obtained 
from 


V _ Z°{T)Z?,{T F ) 
ZZ(T V )Z$(-U) 


(44b) 


The term Z£ is the vibrational partition function for 
the dissociating species and kj r is the dissociation 
rate constant that would exist under conditions of 
thermal equilibrium. The temperature Tp is defined 
as 


J_ _ Jl _ I _ 1 

7 > " Y v ~ T U 


(44c) 


where the quantity — U may be considered as the 
vibrational temperature at which the molecules are 
formed by recombination. The negative value relates 
to the fact that, on the basis of an exponential 
distribution, more molecules are formed in upper 
vibrational levels than in lower levels. Marrone and 
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Treanor note that a value of U ~ Ef^ r /3k gives good 
comparisons between experiment and computation 
for dissociation lag times behind a shock. (See 
tables III and IV for values of activation energy E/ r .) 

Park model (ref. 32). In this model Park assumes 
that certain classes of reactions can be described by 
a single rate-controlling temperature which is an ap- 
propriate average of the local translational, vibra- 
tional, and electronic temperatures. He suggests 
the use of a temperature (weighted heavily with 
the vibrational temperature due to the preferential 
dissociation concept) defined by 

T d = (TTy) 1 / 2 (45) 

to characterize dissociative reactions. This relation 
is empirical and has produced good comparisons 
with some experimental data for radiative energy 
flux (ref. 32), but the model cannot be justified on 
this basis alone. However, it does reproduce cor- 
rect phenomenological trends and is simpler to ap- 
ply than the corrections given by equations (44). A 
more recent investigation (ref. 59) found that a mi- 
nor variation of equation (45) in which T ^ = T ,7 T^ 3 
gave results for reaction rate coefficients that were 
within a factor of 3 of those calculated on the ba- 
sis of the theory of Schwartz, Slawsky, and Herzfeld 
(SSH theory, ref. 60). It is also assumed that the 
rate-controlling temperature for electron impact ion- 
ization is Ty ( T e for a three-temperature model) be- 
cause the free-electron translational energy and elec- 
tronic excitation energy characterize these reactions. 
All other reactions are characterized by the heavy- 
particle translational-rotational temperature T. 

The forward and backward reaction rate coeffi- 
cients can now be expressed by 


equilibrium constant of the form 
K c , r = exp(B r 1 +B r 2 lnZ+B r 3 Z+B r 4 Z 2 +B r 5 Z 3 ) (47) 
where 

Z = 10000/T (48) 

and the constants B J are presented in table III. 

Reaction Sets and Reaction Rate Coefficients 

Two sets of chemical reactions and reaction rate 
coefficients have been employed within the context 
of a two-temperature environment. There is no 
overwhelming evidence at this time to prefer one over 
the other, so both are presented. 

Park’s proposed set of chemical reactions and re- 
action rate coefficients for his two-temperature model 
(ref. 32) are presented in table III as outlined in the 
previous section. He has provided a set of guide- 
lines for defining the rate-controlling temperature in 
different types of reactions. These same guidelines 
have been applied to another set of chemical reac- 
tions and reaction rate coefficients proposed by Dunn 
and Kang (ref. 25). The list of reactions and asso- 
ciated parameters for Dunn and Kang’s chemical ki- 
netic model are presented in table IV. This model 
was originally presented in the context of a single 
temperature, but Park’s guidelines for defining the 
rate-controlling temperature in dissociative and elec- 
tron impact ionization reactions have been employed. 
Dunn and Kang defined the backward rate coefficient 
directly in the form 

h,r = C b , r T nb ’ r exp {-E b ,/kT) (49) 

The parameters needed to define equation (49) are 
included in table IV. 


k f,r = c f, r Tq f ' r exp(-Ef >r /kT q ) (46a) 


k b,r = 


fc/,r (T) 

K c , r 


(46b) 


where K c , r is the equilibrium constant for the r 
reaction and the preexponential parameters Cf^ r and 
Uf r and the activation energy Ej r divided by the 
Boltzmann constant k are presented in table III. The 
reactions, from which /3 s?r and a SyT can be deduced, 
are also presented in table III. The term T q is a 
dummy variable for the rate-controlling temperature. 
The equilibrium constant can be determined from 
the activation energy of the forward reaction and 
the partition functions of the reactants and products 
(ref. 52). Park (ref. 61) employed a curve fit for the 


Vibrational Energy Reactive Source Terms 

The variable which appears in equations (3) 
and (16), represents the vibrational energy per unit 
mass of the diatomic molecules, which are created 
or destroyed at rate w s . If one assumes preferential 
dissociation and recombination of molecules in the 
higher vibrational states (i.e., a molecule is more 
likely to dissociate if it is in a higher vibrational state 
and atoms that recombine are more likely to create 
molecules in a higher vibrational state), then D s 
should be larger than the average vibrational energy 
e v , s or ey s of the system. The value of D s could be 
taken as some fraction of the dissociation energy of 
the molecule, D s . Thus, 


D s — c\D s (50a) 
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where c\ is a constant less than 1 and values for D s 
are provided in table I. Park (ref. 62) has recently 
suggested a similarly motivated definition, 

D s — D s — kT (50b) 

which assumes that the vibrational energy removed 
by dissociation of one molecule comes from an energy 
level which differs from the dissociation energy by the 
average translational energy. A value of c\ = 0.8 has 
been used in equation (50a) for some test calcula- 
tions. (Recent theoretical work of Sharma, Huo, and 
Park (ref. 59) indicates a value of c\ = 0.3 may be 
more appropriate.) In practice, the application of ei- 
ther equation (50a) (with ci = 0.8) or equation (50b) 
tends to lower substantially the vibrational tempera- 
ture behind the shock where dissociation occurs and 
to raise substantially the vibrational temperature in 
the boundary layer where recombination occurs when 
compared with a calculation using D s = ey )S . The 
substantial changes are caused by the large difference 
between D s and ey )S in the test calculation when D s 
is a factor of 10 or more greater than ey >5 across the 
shock layer. In light of this result, a more moderate 
approximation (though still empirical) is to define 

D s - c 2 ^V,s ( 51 ) 

where £2 is a constant greater than 1 for the prefer- 
ential dissociation model and equal to 1 for a non- 
preferential model. This approximation implies that 
the preferentially dissociated molecules are generally 
in higher vibrational states than the average molecule 
but need not be near the dissociation energy before 
a collision in order to dissociate. 

Electronic Energy Reactive Source Terms 

The term he, sis in equations (4) and (16) 

accounts for the rate of electron energy loss when a 
free electron strikes a neutral particle and frees an- 
other electron (ionizes the particle), with a resulting 
loss in electron translational energy. The subscripts 
6 < s < 10 account for the five ionized species which 
potentially should be considered, n e , 5 is the molar 
rate of ionization producing these species, and I s is 
the first ionization energy of the species per mole 
and is presented in table I. The only reactions con- 
sidered in either the Park model or the Dunn and 
Kang model involving electron impact ionization are 

O + e“ « — > 0 + + e~ + e" 

N + e" 4 — > N + + e~ + e~ 

The molar rate of ionization is simply the forward 


reaction rate where subscript r refers to the 

reactions listed above. Thus, 

™e,s =: -R/,r 

where for s = 6 , r = 21 for the Park model (table III) 
and r — 22 for the Dunn and Kang model (table IV), 
and for s — 7, r — 20 for the Park model and r = 21 
for the Dunn and Kang model. 

Note that this model assumes that all of the en- 
ergy required to ionize the species comes from elec- 
tron translational energy and the ionization energy 
is taken from the ground state. This probably over- 
estimates the electronic energy loss rate due to elec- 
tron impact ionization. In the test problems consid- 
ered herein, this term is a small contributor to the 
overall energy balance in the two-temperature model. 

Nonpreferential Dissociation and 

Recombination 

The expressions for the reaction rate coefficients 
for dissociative reactions (eqs. (45) and (46)) and the 
definition of D s used in the vibrational energy reac- 
tive source term (eqs. (50) and (51)) are based on 
the concept, explained previously, that molecules in 
the higher vibrational modes are preferentially disso- 
ciated. However, Jaffe (ref. 63) used collision theory 
to evaluate the reaction rate coefficients for disso- 
ciative reactions in a multitemperature environment 
and predicted a much weaker dependence of on 
T v than that predicted by Park’s model. It may be 
argued that in situations when the free-stream ki- 
netic energy per unit mass is much larger than the 
dissociation energies of the molecules, a vibrational 
ladder-climbing process is not required to dissociate 
the molecules. The methods of statistical mechanics 
give no preferential weighting to any particular en- 
ergy mode in the determination of the total energy 
available in a collision and whether a collision will 
result in an elastic, inelastic, or reactive encounter. 
On this basis, the value of 62 in equation (51) should 
equal 1.0 and the value of T q used in the Arrhenius 
term (exponential term) of equation (45) should be 
weighted by the energy available in all modes. For 
example, 


m ~ e tT + eyTy 

J-a — ~ . - 

(53) 

etr + ey 

Psl e sl,t + Ps2 e s2,t 

e tr — 

(54a) 

Ps 1 “T Ps2 

Psl^sl.V + Ps2^s2,V 

e v — 

Ps 1 + Ps2 

(54b) 
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and si and s2 are the indices of the reactants, e s ,t 
is the translational and rotational energy per unit 
mass of species s, and e s y is the vibrational and 
electronic energy per unit mass of species s. (In or- 
der to accommodate proper weighting at low tem- 
peratures (i.e., T ~ 300) it would be necessary to 
evaluate e Sj t and e s y with respect to a reference 
temperature of 0 K.) This specification has not yet 
been tested; however, sample calculations with either 
Park’s geometric mean temperature (eq. (45)) or the 
translational temperature for dissociative rate con- 
stants give species energy distributions behind the 
shock that cause the value of T q in equation (53) to 
be heavily weighted toward the translational temper- 
ature T. 

The models developed from preferential and non- 
preferential dissociation assumptions are empirical 
and are indicative of the uncertainty in the details of 
the kinetics. Further discussion of these points can 
be found in references 58, 60, and 64. Comprehensive 
quantum-mechanical theoretical studies and non- 
obtrusive laser diagnostic experimental studies are 
required to refine these models. 

Relaxation Processes 


Vibrational-Translational Energy Relaxation 

Millikan and White (ref. 65) present semiempir- 
ical correlations between observed vibrational re- 
laxation times over a temperature range of 300 to 
8000 K and the relevant molecular constants. These 
correlations permit an estimation of r ^ W , the vibra- 
tional relaxation time for species s due to inelastic 
collisions. The correlation is expressed as 

^ rij exp [a. (r' 1/3 - 0/015 /x^ 4 ) - 18.42] 

MW t=l 

pi's To 

E n 3 
j= 1 

(55) 

where p s j is the reduced molecular weight of the 
colliding species s and j,rij is the number density 
of species j, and p is in atmospheres. Values of A s 
for different molecules s are given in table I. (More 
recent correlations, valid over a temperature range of 
300 to 9000 K, are presented in ref. 66.) 

For temperatures above 8000 K, Park (ref. 29) 
suggests an expression for the vibrational relaxation 
time of the form 


rf = (a s c s n s ) 1 (56) 


where c s , the average molecular velocity of molecule 
s, is expressed by 

c 5 - (SkT/irnis) 1 ^ 2 (57) 

n s is the number density of molecule s, and a s is 
the effective cross section for vibrational relaxation. 
(Park noted that eq. (55) yields cross sections for 
vibrational relaxation that are far too large at tem- 
peratures above 8000 K.) The effective cross section 
is assumed to be an order of magnitude smaller^than 
the elastic cross section. It is set equal to 10“ cm 
in the present model. The blending of the two 
relations is accomplished by defining 

<r s >=r^ + rf (58) 

Park has shown (ref. 29) that the variation of t s 
given by equation (58) agrees better with the avail- 
able data for experimental O 2 vibrational relaxation 
time over a temperature range of 5000 to 8000 K 
than does the correlation of Millikan and White given 
by equation (55) . More data are needed to better 
model the relaxation processes at high temperatures 
(T > 8000 K). 

The vibrational energy relaxation terms in equa- 
tions (3) and (16) can be further simplified by making 
the following approximations: 

e* v , s -e v , s *C s v>v (T- T v ) 

E PsC'VyV ^ pC*V,V 

, < T S > ~ Ty 

s= mol. 


(59) 

(60) 


where 


Tv 


E Ps/ M S <T S > 

s— mol. 

E Ps /Ms 

mol. 


(61) 


These relations reduce the number of species- 
dependent parameters which must be evaluated and 
carried along in the calculation. They also permit 
the vibrational relaxation terms to be evaluated as a 
function of a single relaxation coefficient, given by 
equation (60), times the difference in the transla- 
tional and vibrational temperatures. These approxi- 
mations are believed to be consistent with the current 
model, which specifies a single vibrational tempera- 
ture for all molecules. The temperature difference, 
which drives the relaxation process, is treated im- 
plicitly in the numerical algorithm. Test calculations 
show that the vibrational relaxation coefficient can 
be treated explicitly and time lagged, even at large 
Courant numbers, and still achieve stable, convergent 
solutions. 
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Note that the simple Landau- Teller equation for 
expressing vibrational relaxation rates used in equa- 
tions (3) and (16) assumes that the relaxation rate 
varies linearly with the difference in the vibrational 
energies e* — e v . At high temperatures, vibrational 
relaxation obeys a diffusion-like equation with re- 
spect to vibrational energy levels (ref. 32), and a cor- 
rection is needed for the vibrational relaxation time 
given by equation (61). Park (ref. 32) suggests an 
appropriate modification of the form 


1 = ± T sh - T v s 1 

f v Tv T sh - T v sb 


where 


a = 3 . 5 exp ( — 5000 /T sh ) 


(62) 

(63) 


and ?sh and sh are respectively the translational- 
rotational and vibrational temperatures at the point 
on the shock where the relaxation process was ini- 
tiated. (A recent paper (ref. 67) provides species- 
dependent characteristic temperatures for the 
argument of the exponential in eq. (63).) In one- 
dimensional flow, the determination of the post- 
shock temperatures is trivial. In two- or three- 
dimensional flow, a rigorous treatment of the model 
calls for tracing the streamline back to the point of 
origin at the shock. However, within the context 
of the total approximate nature of this model, it is 
more appropriate to choose some average post-shock 
temperature for multidimensional flows. The impor- 
tance of this correction, as seen from the bridging 
formula of equation (62), depends on the magnitude 
of the average post-shock translational-rotational 
temperature. 

Electronic-Translational Energy Relaxation 


Term 6 of equation (4) and term 7 of equa- 
tion (16) were derived by Appleton and Bray (ref. 28) 
to model the energy exchange for elastic collisions 
between electrons and atoms and between electrons 
and ions. The frictional heating of electrons by 
heavy particles due to differences between electron 
and heavy-particle velocities is ignored because of 
the assumption of ambipolar diffusion. Appleton and 
Bray’s model was for plasmas which did not contain 
molecules with internal degrees of freedom. Their 
expressions for u es , the effective collision frequency 
of an electron with species s, are presented below 
for collision partners being either neutral or charged. 
For Coulomb collisions between electrons and ions, 
the expression is 


Ves 


8 

3 


/M 1 / 2 4 1 

\m e ) (2feT e ) 3 /2 


In 



(64) 


where n e is the number density of electrons, n s is the 
number density of species s, m e is the electron mass, 
e is the magnitude of the electronic charge equal to 
4.80298 X 10 -10 esu, and k is Boltzmann’s constant. 
For collisions between electrons and neutrals, the 
effective collision frequency is expressed as 


L'es — TlsCTes 


/s fer e \ 1/2 

V 7T m € ) 


(65) 


where the effective electron-neutral energy exchange 
cross section is defined by a curve fit of the form 


(Tes — a s + b s T e + c s Tg ( 66 ) 

The constants for equation (66) are presented in ta- 
ble V. This curve fit was generated from effective col- 
lision cross-section data at 5000, 10 000, and 15 000 K 
found in reference 68. 

Vibrational-Electronic Energy Relaxation 

Term 6 of equation (3) and term 8 of equa- 
tion (4) in the three-temperature model call for 
an approximation to the effective relaxation time 
for vibrational-electronic energy accommodation 
< T es >. Recent data of this type are discussed by 
Lee (ref. 69) and Huo et al. (ref. 70). Much of these 
data are in tabular form, and no attempt was made 
to curve fit the data because of the emphasis on the 
two- temperature model. 

Transport Properties 

The derivation of the transport properties in this 
section closely follows the example of Lee (ref. 2). 
The approach is based on an extension of Yos’ for- 
mula (ref. 71) to the multitemperature gas mixture. 
The collision integrals for heavy particles are based 
on the heavy-particle translational temperature T. 
The collision integrals for electrons with any other 
partner are based on the electron temperatures T e or 
T v . The collision integrals are evaluated as curve fits 
to the tabular data generated in reference 68. The 

curve fits assume a linear variation of log 10 

and log 10 (ttQ ^’ 2) ) with ln(T), as defined by the 

data for these quantities at 2000 and 4000 K pre- 
sented in table VI. Therefore, 

log i 0 (-n^ fc) ) 

= log 10 (nni k r ' k) ^ (2000) 

loglO (4000) - log 10 (*-n ir’ k) ^ (2000) 

+ ln(4000) - ln(2000) 

X [ln(T) - ln(2000)] 
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This temperature range is used to give the best over- 
all agreement within the boundary layer for typical 
flow-field simulations. In equation (67), T e or Ty is 
used for T in collisions involving electrons. The units 
of n sr are square centimeters. The data used in ta- 
ble VI assume an electron pressure of 0.001 atm. In 
general, the collision integrals should be corrected for 
the effects of varying electron pressure as discussed 
by Armaly and Sutton (refs. 72 and 73). The fol- 
lowing formula of reference 68 may be employed to 
correct the collision integrals tabulated in table VI 
for a different electron pressure: 

*tffi fc) (io- 3 ) 

ni nir’ fe) (Pe) 

ln[20.9(r/1000) 4 + 152(r/1000) 8 / 3 ] /ggx 
~ ln[0.0209(T 4 /10 12 p e ) + 1.52(T 4 /10 12 p e ) 2 / 3 ] 


where p e is the electron pressure in atmospheres. 

Two modified collision integrals which are used 
extensively in subsequent evaluations of transport 
properties are defined below. 


A «(T) = | 

A^(T) = f 


2 MsMr 

nRT{Ms + Mr) . 

2 MsMr _ 

7rlT(M s + M r )\ 


I 1/2 




1/2 


7rn£ ,2) 


(69) 

(70) 


The molar concentration of species s, q s , is also used 
extensively in subsequent calculations. It is defined 

as 

Ps 


Is = 


pM s 


(71) 


The mixture viscosity // can now be expressed as 


10 


-E 

S = 1 



The translational energy thermal conductivity of 
heavy particles rjt is expressed as 


10 


Is 


m ~ 4"^ 10 ... 

» =1 r'yrA^(T) + 3.54'y e A^ ; (T e ) 

r — 1 


,(2)l 


(73) 


where a sr is defined by 


a sr 


[1 — (m 3 /m r )][0.45 — 2.54(m a /m r )] 
~ 1 + [1 + (m 8 /m r ) l 2 


The rotational modes of molecules are assumed 
to be fully excited and the rotational energy thermal 
conductivity r] r is expressed as 


Vr 


* E 

s-mol. 


Is 

2'rrA^ ) (T)+ 7e Aie ) (Te) 
r = 1 


(75) 


The frozen thermal conductivity of the mixture for 
translational and rotational energy of heavy particles 
r] is now given by 


V = m + r\ T 


(76) 


The vibrational thermal conductivity r] v equals the 
rotational thermal conductivity T] r . Thus, 


Vv = Vr 


(77) 


The electron thermal conductivity r) e follows the 
form of equation (73) and is given by 


15 k 'll 

4 E 1.45q r A^(T e ) 
r= 1 


(78) 


The binary diffusion coefficient for a pair of heavy 
particles sr is defined as 


D 


sr 


kT 

pAil\T) 


(79) 


The binary diffusion coefficient between electrons 
and heavy particles is expressed as 


D 


er 


kTe 

P A { Jr\T e ) 


(80) 


The effective diffusion coefficient of species s in the 
mixture, D s , used in the governing equations can now 
be evaluated by 


7, 2 M s ( 1 - Msis) 
E (lr/D sr ) 

r— 1 
r^s 


(81) 


where 7 1 is defined by 


7 1 - Is (82) 

S— 1 

The diffusion of ions and electrons is linked be- 
cause of the induced electric field which occurs in 
the presence of an electron pressure gradient. In a 
partially ionized gas with zero electric current, this 
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effect is modeled with the ambipolar diffusion coeffi- 
cient D? on , where 

D ion = 2 Aon (83) 

Each ion in a multicomponent gas mixture is assumed 
to diffuse as if it were the only ionic species in the 
mixture. Therefore, the effective diffusion coefficient 
of ions is set equal to the ambipolar diffusion coeffi- 
cient as defined in equation (83). Within the context 
of the ambipolar diffusion approximation, the effec- 
tive diffusion coefficient of electrons D e is obtained 


by equating the diffusion velocity of electrons with 
the diffusion velocity of ions. This specification leads 
to the following relation for D e : 

10 

E d- s1s 

D e = m e ~-~- (84) 

E msls 

s=6 

The relations in this section completely define 
all the transport properties used in the governing 
equations. 


Upwind Formulation of the Flux Vector 


Upwind, or total variation diminishing, numerical formulations of the governing conservation laws have 
been shown to be robust with regard to their capabilities to simulate hypersonic flows with strong shock 
waves and generally complex wave interactions (refs. 35, 36, 74, and 75). These formulations usually require 
a factorization of the Jacobian of the inviscid flux vector involving the right and left eigenvectors and the 
eigenvalues. In particular, if f is the inviscid flux vector and q is the vector of conserved variables, then the 
Jacobian of f with respect to q is given by 


dq ^ LAR (85) 

where A is a diagonal matrix composed of the eigenvalues of A, L is a matrix of column eigenvectors, R is a 
matrix of row eigenvectors, and LR equals the identity matrix I. 

It is convenient to formulate the numerical solution of the governing equations within the context of a 
finite- volume scheme. The finite-volume schemes work with the integral forms of the conservation laws and 
set up approximations to flux across cell walls defined by the distribution of neighboring grid points. Once 
the cell wall is defined, it is a trivial procedure to define unit vectors which are normal and tangent to the cell 
wall. The extension to a finite- difference scheme is obtained by noting the relationship between the ratio of 
cell wall areas to cell volumes in the finite-volume formulations and the metric coefficients in finite-difference 
formulations (ref. 76). 

The vector of conserved variables in the two-temperature model defined by equations (1), (2), (5), and (16) 
within the context of a finite-volume approximation is presented below. 


Ps 

pu 



pE 


Ipey J 

The inviscid flux vector for the two-temperature model is written as 


( 86 ) 


pU c s 

pUu + pn x 
pU v + pn y 
pU w + pn z 
pUH 
pU e v 


( 87 ) 
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where c s is the mass fraction of species s, defined by 


C S = ^ 

p 


(88) 


n x ,ny, and n z are the x,y, and z components of a unit vector normal to a computational cell face, and U is 
the normal component of velocity through the cell face, defined by 


U = un x + vn y + wn z 


(89) 


It is useful in the evaluation of the eigenvectors to employ two unit vectors, 1 and m, such that n, 1, and m 
are mutually orthogonal (i.e., — l^rui = 0). The velocity components in the 1 and m directions, 

tangent to the cell face, are then defined by 


V = ul x + vly + wl z 
W = um x + vrriy + wm z 

The Jacobian of f with respect to q is expressed as 



U (6 sr c s ) 

c s n x 

C S Uy 

c s n z 

0 

0 


i r n x -Uu 

—f3un x + un x + U 

—f3vn x + un y 

—f3wn x + un z 

P n x 

4>fix 

A - 

IrKy — Uv 

~(3un y 4- vn x 

—f3vn y + vn y + U 

—ftwriy + vn z 

pn y 

(f)Tly 


7 r n z - Uw 

—j3un z + vm x 

—(3vn z + wn y 

~/3wn z + wn z + U 

P n z 

<t>n x 


7 r U - UH 

—(3uU + Hn x 

—(3vU + Hn y 

—fjwU 4 - Hn z 

I3U + U 

4>u 


. —Uey 

eyn x 

eyriy 

e v n 2 

0 

u 

The similarity transformation matrices R and L are defined as 





R = 


L = 


' o?8 sr - Cs^r /3uc s 

(3vc s 

f3wc s 

-0c s 

1 

■s- 

c* 

Co 

1 


-V l x 

ly 

lz 

0 

0 


-W m x 

m 

'y 

m z 

0 

0 


7 r ~Ua an x — (3u 

an y 

- f3v 

an z — (3w 

P 

</> 


7 r 4* Ua — an x — f3u 

an y 

— f3v 

—an z — (3w 

P 



—ey^fr fiuey 

flvey 

f3wey 

-Pev 

a? — 4>ey _ 


6 sr /a 2 

0 

0 

c s /2a 2 

c 

■s/2 a 2 

0 ‘ 

u/a 2 

lx 

m x 

(u 4- an x )/2a 2 

[u- 

an x )/2a? 

0 

v/a 2 

ly 

m y 

{v 4- an y ) /2a 2 

{v- 

an y ) /2a? 

0 

w/a 2 

lz 

m z 

(w + an z )/2a 2 

(w — 

an z )/2a? 

0 

[/3 {u 2 +v 2 + w 2 ) — 7 r ]//3a 2 

V 

W 

(H + all) /2a 2 

(H- 

- aU)/2a? 

-4>/f3a 2 

0 

0 

0 

ey /2a 2 

ey /2a? 

1/a? 


(90) 

(91) 


(92) 


(93) 


(94) 


21 



The diagonal matrix of eigenvalues of A is defined by 


'U 0 0 0 0 0 | 

0 U 0 0 0 0 

0 0 U 0 0 0 

A = (95) 

0 0 0 U + a 0 0 

0 0 0 0 U — a 0 

.0 0 0 0 0 U . 

In the matrices defined above, the first row and column correspond to the 11-species continuity equations. 
Subscript s refers to row s and species s and subscript r refers to column r and species r, where both s and 
r vary from 1 to 11 in the present model. The variables f3, <j>, and qy are related to the partial derivatives of 
pressure with respect to q, and a is the frozen speed of sound. These quantities are derived below. 

The differential form of equations (8) and (9) for pressure can be written as 



10 

dp = RTY, 
8=1 

dps | -Jg Y P s | dp e J R dTy p e 

M s ^ M s M e M e 

8=1 

(96) 

We need to express dp as a function of dq. This relation can be established by first expressing dT and dTy 
as a function of de and dey with the equations presented in the Thermodynamic Relations section. The two- 
temperature corollary to equation (32), written with respect to energy e as opposed to enthalpy h , is expressed 
as 



f T 

es — Cy tr dT' + ey )S + e s>0 

■IT ref 

(97) 

and 


rTy 

ev,s = / C° v y dT' 

(98) 

where 


Si 

II 

+ 

•s 

(99a) 



/mS /OS | /^yS 

^ VfV ^VjV ' 

(99b) 

It is convenient to define the specific heat capacity at constant volume for a free electron as C% y and to set 
Cy tr = 0 because the energy of an electron is a function only of Ty. The differential forms of equations (97) 
and (98) are given by 

de s = C s V)tr dT + C s vy dTy (100) 

and 


dev,s = dT v 

(101) 

Recall that 


11 

e = c « e s 

S=1 

(102) 

and 


11 

e V °s e V,s 

8 = 1 

(103) 
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so the differential expressions for de and dey can be written 


11 


11 


11 


de= ^2 dc s e s + E Cs des = E dc * e s + E c s{Cy ttr dT + C*y dT v ) 


S— 1 


11 


S=1 


11 


s=l 


11 

e 

5—1 


11 


de V - E dc * e v,s + E Cs de V,s = E dCs e V,s + E CsC v,V dT V 


5=1 5=1 S=1 

Solve for dT and dTy with equations (104) and (105) to obtain 


11 

E 

S = 1 


(104) 


(105) 


li 


dT = 


de — dey — ^2 dc s {e s — ey s ) 


C v , } 


tr 


11 


dey — Y2 dc s ey £ 


dTy == 


5=1 


' VyV 


(106) 


(107) 


where equations (39) and (88) have been used to obtain the heat capacities of the mixture. The differential 
form of energy de, vibrational energy dey, and mass fraction dc s can be written with respect to the elements 
of dq from equations (13), (86), and (88) as follows: 


de = 
dey = 


_ dpE — E dp — (u dpu + v dpv + w dpw ) + ( u 2 + v 2 + vr)dp 

P 

dpey — ey dp 


dc. s — 


dps Cg dp 


(108) 

(109) 

( 110 ) 


Note that dp is easily expressed in terms of the elements of dq by 

11 

«*/> = E d f jg 

s= 1 


(111) 


Substitute equations (106) to (111) into equation (96) and combine terms to obtain 


dp — f3(dpE — u dpu — v dpv — w dpw) + </> dpey + 2s dp s 


where 


and 


Q — dP _ ^ y' Pr 

dpE pC V £ r M r 

4 , = _Sp_ = JE _^_ /3 

dpey pC v y M e 


Is - 


dp 

dp s 


RTq U 2 +V 2 +W 2 

M s 2 


- /3e s - (j>e v , s 


(112) 


(113) 

(114) 


(115) 
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In equation ( 115 ), T q = Ty when 5 is an electron; otherwise, T q = T . The frozen speed of sound a can now be 
evaluated with equations (113) to (115): 

a 2 = c 5 7 s + P[H - (u 2 + v 2 + w 2 )] + (j>e v = (1 + /?)- ( 116 ) 

s=l 9 

This definition of a 2 comes from the evaluation of the eigenvalues of A. 


Results and Discussion 

The task of validating the collection of physical 
models assembled here with regard to their predic- 
tive capabilities for hypersonic flows in chemical and 
thermal nonequilibrium is beyond the scope of this 
paper. Elements of the models have been validated 
to a limited extent in the original sources. How- 
ever, validation of the ensemble in realistic condi- 
tions is currently limited because of the difficulty in 
obtaining the experimental data. It is expected that, 
as validation studies proceed using some data now 
available (refs. 77 and 78) and data which may be 
available in the near future (ref. 79), the models will 
evolve from the analytical forms or parametric curve 
fits presented herein. The present paper serves as a 
single source benchmark for these studies and pro- 
vides guidance for implementation of the models in 
detailed computer codes. 

Even without experimental data, it is still possible 
to compare results obtained from options presented 
herein to demonstrate the impact of uncertainties in 
various elements of the model. To this end, some pre- 
dictions have been prepared using LAURA (refs. 35 
and 36) with the strong implicit coupling (ref. 15) 
provided through equations (85) and (93) to (95). 
Comparisons with predictions made by a Direct- 
Simulation Monte Carlo (DSMC) algorithm (ref. 37), 
a kinetic-theory-based particle simulation approach 
to hypersonic, rarefied- flow analysis, provide addi- 
tional opportunities to evaluate present capabilities. 
As with any continuum-based approximation scheme, 
the DSMC approach is also subject to physical mod- 
eling errors, particularly with regard to the descrip- 
tion of real, reacting gases. However, the DSMC 
method is better at describing low-density flows and 
it does not require any a priori assumptions concern- 
ing the evaluation of dissipative phenomena as re- 
quired in the Navier-Stokes approximation. These 
characteristics make DSMC a valuable benchmark 
for evaluating continuum-based Navier-Stokes flow- 
field solutions in transitional flow regimes. 

Chemical Kinetic Model Studies 

Profiles of mixture density, pressure, species num- 
ber densities, and temperatures across the shock 


layer near the stagnation streamline over an axisym- 
metric approximation to the Aeroassist Flight Ex- 
periment (AFE) (refs. 37 and 79) are presented in 
figures 1(a) to l(o). The free-stream conditions are 
free-stream velocity of 8917 m/s, free-stream den- 
sity of 0.0000272 kg/m 3 , and free-stream tempera- 
ture of 197 K. These conditions correspond to an 
altitude of 78 km. Boundary conditions correspond- 
ing to a noncatalytic, no-slip cold wall under zero 
normal pressure gradient have been imposed. The 
grid is defined by 64 cells stretching from the body, 
across the captured shock, and into the free stream 
and 39 cells from the axis of symmetry to the circular 
shoulder. (Grid structure is discussed in more detail 
at the end of this section.) Comparisons are made 
between the predictions obtained with the chemical 
kinetic model of Dunn and Kang (D & K, ref. 25) 
(cases I and II) and the chemical kinetic model of 
Park (ref. 32) (case III). Equation (45) is used to de- 
fine the rate-controlling temperatures for dissociation 
in both chemical kinetic models (cases I and III). In 
addition, a specification of — T is tested within 
the context of the Dunn and Kang chemical kinetic 
model (case II). The rate-controlling temperatures 
for recombination are set equal to the rate-controlling 
temperatures for dissociation. In general, the trans- 
lational temperature should be used for such reac- 
tions; however, the translational temperature is very 
nearly equal to the vibrational temperature where re- 
combination is significant in all test cases presented. 
Equation (51) is used to model the vibrational en- 
ergy lost through dissociative reactions, with c<i = 1. 
The corrections defined by equations (62) and (63) 
are not applied to the calculation of the vibrational- 
translational relaxation time. 

The mixture density and pressure are relatively 
insensitive to variations in the chemical kinetic mod- 
els tested, as shown in figures 1(a) and 1(b). Only 
a slight variation in shock standoff distance is ob- 
served. Chemical nonequilibrium effects are shown 
in the slight overshoot of both neutral and ionized 
molecules behind the captured shock (figs. 1(e) to 
1(g) and l(j) to 1(1)). Specification of the geometric 
average temperature as the rate-controlling temper- 
ature inhibits dissociation, as expected, when com- 
pared with use of the heavy- particle translational 


24 



temperature (figs. 1(e) to 1(g)). There is significant 
formation of O 2 and NO occurring in the boundary 
layer at this condition, even with the noncatalytic 
wall boundary condition, as shown in figures 1(f) and 
1(g). The greatest differences among the three chem- 
ical kinetic models are observed in their predictions 
of ionized-species profiles. Electron number density 
appears to arise from atomic oxygen over most of the 
shock layer (figs. l(i) and l(m)), with significant con- 
tributions from atomic nitrogen immediately behind 
the captured shock. Ionized molecular oxygen pro- 
duced behind the shock is of the order of 0.01 percent 
of the total number density, but it quickly falls off as 
the molecule dissociates and accommodates to local 
conditions (fig. l(k)). Species mole fraction profiles 
across the shock layer from case I for neutrals (fig. 2) 
and for ions (fig. 3) are presented in order to highlight 
the relative concentrations of the constituent species. 

The chemical kinetic models of Dunn and Kang 
predict a significantly different profile for the deion- 
ization of N + and 0 + in the boundary layer than the 
model of Park (figs. 1(h) and l(i)). The exact cause 
of these variations between the two chemical kinetic 
models is unknown at present, but both the reaction 
sets describing charge exchange and ionization and 
the reaction rate coefficients in these sets are very 
different. 

Thermal nonequilibrium is evident in figures l(n) 
and l(o). The spike in translational temperature 
behind the captured shock is indicative of the de- 
layed dissociation due to chemical nonequilibrium, 
and the corresponding low vibrational temperature 
shows the significant thermal nonequilibrium in this 
region. The specification of translational tempera- 
ture as the rate-controlling temperature for dissoci- 
ation yields lower peak translational temperatures, 
enhanced dissociation, and a more gradual accom- 
modation of the vibrational temperature with the 
translational temperature. 

Figure 4 presents the DSMC results, in addition 
to results for the other models, for species mole frac- 
tion and temperature profiles across the shock layer. 
Shock standoff distances for the continuum-based 
LAURA and the noncontinuum-based DSMC algo- 
rithm are approximately equal, as judged by the lo- 
cation of the peak in translational temperature. The 
DSMC shock thickness exceeds the LAURA predic- 
tion by about a factor of 3 as judged by the high gra- 
dient region in mole fraction of atomic nitrogen and 
oxygen spanning from y = 0.001 to post-shock levels. 
The DSMC methods generally predict thicker shocks 
than continuum-based methods and are expected to 
be more accurate than continuum-based predictions 
in this regard at low densities because of their ability 
to better simulate the random motions of particles 


across the shock front. At present, continuum-based 
methods compute this dissipation of mass, momen- 
tum, and energy due to random thermal motions of 
particles as a linear function of gradients in the flow 
field. They have no mechanism to recognize the sig- 
nificant change in mean free path across the shock 
front and generally utilize computational cells that 
are less than a mean free path in length ahead of the 
shock for hypersonic flows at 80 km and above. 

The DSMC results show a separate rotational 
temperature in addition to the translational and vi- 
brational temperatures (figs. 4(e) to 4(g)). These 
kinds of data arise from the modeling of diatomic 
molecules in the system and the monitoring of trans- 
lational, rotational, and vibrational energies of such 
particles. The continuum prediction of peak heavy- 
particle temperature T, which includes both trans- 
lational and rotational energies, falls between the 
peaks of the DSMC predictions for translational- 
rotational temperatures in case I and case III, as 
should be expected. The Park model (case III) gives 
the best overall agreement with temperature dis- 
tribution across the shock layer as compared with 
DSMC predictions. 

All three continuum chemical kinetic models and 
the DSMC predictions are in generally good agree- 
ment with the post-shock levels of species mole frac- 
tion for atomic nitrogen and oxygen (figs. 4(a) and 
4(b)). Post-shock minimums in molecular nitrogen 
mole fraction vary from 0.04 to 0.08, with the DSMC 
showing the greatest concentration (lowest dissocia- 
tion) of molecular nitrogen (fig. 4(c)). The most sig- 
nificant differences are in the predictions of electron 
number densities (fig. 4(d)). Only the DSMC pre- 
diction shows a well-defined peak in electron number 
density in this semilog plot. The sensitivity of this 
profile to details of the chemical kinetic and thermal 
relaxation models, with regard to both profile shape 
and peak profile values, makes electron number den- 
sity measurements in flight an important contribu- 
tor to help resolve unknowns in the present physical 
models. 

All the continuum solutions are generated on 
identical grids with identical numerical parameters. 
Stagnation point heating values for these three cases 
are 205 kW/m 2 for the Dunn and Kang model with 
T d = {' TTv ) 1 / 2 (case I), 200 kW/m 2 for the Dunn and 
Kang model with T&—T (case II), and 184 kW/m 2 
for the Park model with T ^ = (TTi,) 1 / 2 (case III). 
Stagnation point heating for the DSMC calculation 
is 200 kW/m 2 . In light of the many uncertainties 
with the physical models and the weak dependence of 
computed convective heat transfer rates on numerical 
parameters, the excellent agreement between the two 
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solution techniques is believed to be fortuitous. It is 
expected, based on results of references 18 and 80 and 
unpublished calculations using the present method, 
that larger differences between the flow fields pre- 
dicted by continuum Navier-Stokes and noncontin- 
uum DSMC algorithms will occur across the shock 
transition zone at approximately 90 km, but predic- 
tions of convective heating at the surface may agree 
up through 100 km for AOTV applications. 

Grid Refinement Studies 

The breakdown of the continuum-based methods 
in the transitional flow regime mentioned in the pre- 
vious section leads to subtle contradictions with re- 
gard to evaluating solution accuracy with grid refine- 
ment studies. For example, if the goal is to obtain 
an accurate solution of a linear system of partial- 
differential equations and the approximation scheme 
is stable and consistent, then grid refinement stud- 
ies yield computed (difference) solutions that con- 
verge to the exact solution (Lax’s Equivalence The- 
orem, ref. 81). (This approach is used as well in the 
study of nonlinear systems of conservation laws, but 
care must be used in interpreting results because of 
the possibility of multiple, entropy-violating solution 
branches.) However, if the goal is to obtain an accu- 
rate simulation of physical phenomena under the con- 
straint of a continuum-based Navier-Stokes approx- 
imation scheme, then excessive grid refinement may 
be counterproductive. The Navier-Stokes approxi- 
mation fails to resolve accurately high Mach number 
shock structure. The failure arises from the inad- 
equacy of linear functions of velocity and tempera- 
ture gradients to describe correctly shear stresses and 
conduction in this high gradient region. Cell dimen- 
sions which are larger than a mean free path (cell 
Knudsen numbers smaller than 1) tend to give re- 
sults which agree better with the DSMC calculations. 
Such a restriction on cell size must be viewed strictly 
as an empiricism which better mimics the DSMC 
shock structure and, perhaps, better models the dis- 
sipative phenomena across a shock, given the con- 
straints of the Navier-Stokes approximation. These 
concerns become more acute as altitude increases or 
density decreases. For example, preliminary calcula- 
tions from LAURA with shock-transition-zone Knud- 
sen numbers greater than 1 (not presented herein) 
show a much greater difference in shock thickness 
at 90 km, where the free-stream density is a fac- 
tor of 5 smaller than at 78 km. These problems of 
modeling dissipation with continuum-based methods 
in the transitional region between the free molecu- 
lar and continuum flow regimes will need to be ad- 
dressed more rigorously if routine application of these 
analysis tools to AOTV is to be achieved. 


A grid refinement study was implemented which 
was designed to check the effects of truncation error 
on the computed solutions while maintaining a cell 
Knudsen number less than 1, as discussed above. 
Results of this study for test case I are presented 
in figures 5(a) to 5(j). Grids 1 and 2 are made up 
of 64 cells which are exponentially stretched from 
the body to the inflow boundary ahead of the shock. 
Grid 1 has a minimum cell size at the wall equal 
to 2.872 X 1(T 6 m and an average cell size through 
the shock transition zone equal to 1.5 X 10“ 2 m. 
Grid 2 has a minimum cell size at the wall equal 
to 105 X 10 -4 m and an average cell size through 
the shock transition zone equal to 8.0 X 10 -3 m. 
Grid 3 is made up of 128 exponentially stretched 
cells with a minimum cell size at the wall equal to 
5.12 X 10 3 m and an average cell size through the 
shock transition zone equal to 4.0 X 1(P 3 m. The 
ratio of mean free path to cell size (cell Knudsen 
number) based on conditions at the beginning of the 
shock transition zone equals 0.193 for grid 1, 0.363 
for grid 2, and 0.725 for grid 3. None of the grids 
violates the empirical constraint on minimum cell 
size discussed above, although grid 3 is very close to 
the limit. The cell Reynolds numbers 2 {pa A z/p) 
at the body and in the shock transition zone are 
respectively 0.334 and 26.87 for grid 1, 10.8 and 12.2 
for grid 2, and 5.5 and 5.27 for grid 3. Grid 1 
sacrifices resolution at the shock for a very fine 
resolution through the boundary layer. Grids 2 and 
3 utilize the same exponential stretching parameters. 
Grids 1 to 3 have the same lateral cell distribution 
around the body. 3 The actual distribution of mesh 
points across the shock layer is indicated by the 
symbol location in figures 5(a) to 5(d). 

Figure 5(a) shows a sharper, higher peak in the 
translational temperature as the grid is refined across 
the captured shock. This trend is similar to one ob- 
tained in reference 82, which presents a simulation 
of the effects of a relatively coarse grid resolution of 
reacting flow crossing a normal shock through the 
use of artificial viscosity. The thickness of the peak 


2 The cell Reynolds number, as defined here, is proportional 
to the inverse of the cell Knudsen number because the viscos- 
ity p is proportional to the product of density, sound speed, 
and mean free path (paX). The cell Reynolds number is of- 
ten used to assess adequacy of the grid in the bound ary- layer 
flows. The cell Knudsen number is a natural parameter to 
consider when dealing with low-density flows. Both are doc- 
umented here for the readers’ convenience. 

This study was completed sometime after the calculations 
were made for figs. 1 to 4. A different grid distribution 
function, which tended to an average of grids 1 and 2 in fig. 5, 
was used in the earlier work. 
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(thermal and chemical relaxation zone) for grid 3 
is approximately equal to the DSMC result in fig- 
ure 4(e). Further refinement would probably sharpen 
the peak a little more, but it is not clear that cell 
averages taken over a dimension less than the local 
mean free path would be physically meaningful. The 
temperature gradient approaching the wall is largest 
for grid 1. The vibrational-electronic temperature 
distribution is relatively insensitive to grid once the 
shock is crossed (fig. 5(b)). Differences in density 
(fig. 5(c)) relate mostly to the degree of dissociation 
of nitrogen and to the location of the shock transition 
zone. Post-shock-transition-zone pressure levels are 
independent of grid but the sharpness of the shock 
front is better resolved with the finest grid (fig. 5(d)). 
The mole fraction of atomic oxygen across the shock 
layer is insensitive to the grid (fig. 5(e)); however, 
the mole fractions of molecular and atomic nitro- 
gen (fig. 5(f) and 5(g)) show some dependence on 
the shock-transition-zone processing of these species. 
Oxygen is fully dissociated at this condition and so 
is not sensitive to the details of the flow through the 
shock transition zone. Nitrogen is not fully dissoci- 
ated and greater dependence on resolution through 
the shock transition zone is to be expected. Post- 
shock-transition-zone levels of molecular and atomic 
nitrogen are in good agreement for grids 2 and 3; but 
there is a trend in which coarse grids promote recom- 
bination in the boundary layer. There is generally 
good agreement across the entire shock layer for the 
predictions of 0 + from grids 2 and 3 (fig. 5(h)). The 
differences here tend to be less than the differences 
caused by unknowns in the physical models discussed 
previously for figure 4(d). Ionized molecular species 
appear only in the shock transition zone in any sig- 
nificant levels (fig. 5(i)). The rapid production and 
depletion of ionized, molecular species in the shock 
transition zone shows up as a small plateau in the 
free electron distribution (fig. 5(j)). 

Concluding Remarks 

The conservation equations for simulating hyper- 
sonic flows in thermal and chemical nonequilibrium 
and details of the associated physical models have 
been presented. These details include the curve fits 
used for defining thermodynamic properties of the 
11-species air model (N, O, N 2 , 02, NO, N + , 0 + , 
Nj , Oj , NO + , e”), the curve fits for collision cross 
sections, the expressions for transport properties, the 
kinetic models, and the vibrational and electronic en- 
ergy relaxation models. The expressions were for- 
mulated in the context of either a two- or three- 
temperature model. Greater emphasis is placed on 
the two-temperature model, in which it is assumed 
that the translational and rotational energy modes 


are in equilibrium at the translational temperature 
and the vibrational, electronic, and electron transla- 
tional energy modes are in equilibrium at the vibra- 
tional temperature. The eigenvalues and eigenvec- 
tors associated with the Jacobian of the flux vector 
have also been presented in order to accommodate 
the “upwind” based numerical solutions of the com- 
plete equation set. Thermodynamic relations involv- 
ing the partial derivatives of pressure with respect 
to the conserved variables were derived within the 
context of the two- temperature approximation. 

Two chemical kinetic models and two prescrip- 
tions for the rate-controlling temperature of dissocia- 
tion were studied and compared with 
Direct-Simulation Monte Carlo (DSMC) predictions 
for hypersonic flow over an axisymmetric approxi- 
mation to the Aeroassist Flight Experiment vehicle. 
Differences among the models range from a factor 
of 2 for degree of dissociation to a factor of 10 for 
degree of ionization. Park’s chemical kinetic model, 
which uses the most recent available kinetic data, is 
in closest agreement with the DSMC results for tem- 
perature distributions across the shock layer. All the 
predictions show that electron number density is bal- 
anced by the ionized atomic oxygen number density 
over most of the shock layer for the test condition. 
Predictions for neutral atomic species are in gener- 
ally good agreement; however, the strong dependence 
of profile shape and magnitude of charged particles, 
particularly electrons, on variations in the kinetic 
models highlight the importance of obtaining more 
theoretical and experimental data at flight conditions 
for an aeroassisted orbital transfer vehicle (AOTV). 
Such data will be used to validate and improve the 
present chemical kinetic and thermal relaxation mod- 
els so that the simulation of hypersonic flows at high 
altitudes (where chemical and thermal nonequilib- 
rium effects are important) can proceed with greater 
confidence. 

A grid refinement study was implemented to 
check the effects of truncation error on the computed 
solutions for one of the cases discussed above. Re- 
finement of the shock transition zone sharpens and 
raises the peak translational temperature but has lit- 
tle influence on the post-shock-transition-zone trans- 
lational or vibrational temperatures. Oxygen was 
fully dissociated in this test case, so the post-shock- 
transition-zone levels of atomic oxygen are insensi- 
tive to grid refinement. The two grids with the finest 
resolution through the shock transition zone predict 
equivalent levels of atomic and molecular nitrogen 
behind the shock. Nitrogen dissociation is inhib- 
ited by the coarsest grid through the shock transi- 
tion zone. In a like manner, coarse grids tend to 
promote recombination in the boundary layer. Grid 



refinement through the shock transition zone may be 
considered excessive when computational cell sizes 
become smaller than the local mean free path. 

Limitations of the equations and models pre- 
sented herein may be categorized as parametric and 
physical. Parametric limitations arise from uncer- 
tainties in modeling a physical process. For example, 
the equation sets discussed herein have the flexibil- 
ity to model approximately the effects of preferential 
dissociation on the computed flow field. However, 
the magnitude of this effect is not fully understood 
at this time. Other parametric limitations include 
the thermodynamic and collision cross-section curve 
fits, particularly at high temperatures, and the chem- 
ical reaction sets. The constants within these models 
and/or the models themselves are likely to require 
adjustment as better data become available for com- 
parison. Ultimately, the accuracy of the simulation 
is a function of the uncertainty in the values of the 
parameters which define the models. 

Physical limitations arise from intentional simpli- 
fications and assumptions made to model the phys- 
ical system. For example, the use of the Navier- 
Stokes equations implies a linear relation between 
velocity gradients and shear stresses which are not 
valid through strong shocks. At high altitudes and 
velocities, the internal shock structure becomes a sig- 
nificant part of the forebody merged layer flow field in 
which the shock and boundary layers overlap. Gen- 
erally good comparisons between computed results 
from the continuum approach described herein and 
those from the kinetic approach which uses DSMC 
were obtained for a representative AOTV trajectory 
point at an altitude of 78 km (Mach 32). 

Other physical limitations include the two- 
temperature approximation and the assumption of 
ambipolar diffusion. Neither of these approxima- 
tions is expected to place any additional constraints 
on AOTV applications because of the preponder- 
ance of molecular nitrogen in the forebody flow field 
(compared with other molecules) and the relatively 
low ionization levels. At higher velocities, typical 
of Martian return, the two-temperature approxima- 
tion may still prove valid over a significant portion 
of the aeropass when the flow is fully dissociated and 
vibrational energy contributions go to zero. Base 
and near-wake flow-field simulations are also of in- 
terest to AOTV designers because of payload protec- 
tion. Cell Knudsen numbers will exceed 1 in this 
low-density region. The same concerns that limit 
Navier-Stokes approximations across strong shocks 
also apply across the low-density free shear layer. 
Comparisons with DSMC calculations and experi- 
mental data are required to fully understand the 
limitations of the continuum analysis in this region. 


Finally, the present model cannot account for 
plasma-dynamic effects. These are not expected to 
be important in a continuum, forebody flow field, 
where the estimated magnetic pressures are orders of 
magnitude smaller than the fluid pressures in AOTV 
applications. The importance of these effects in the 
base flow region, where electrical conductivity may 
be higher, is not known. 

NASA Langley Research Center 
Hampton, VA 23665-5225 
November 8, 1988 
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Table II. Constants for Curve Fits of Thermodynamic Properties 


(a) Neutrals 


Species 

Range 

(a) 

M 

a 2 

a 3 

a 4 

^5 

^6 

Source 

N 

1 

0.2503071E+01 

— 0.2180018E— 04 

0.5420528E— 07 

— 0.5647560E— 10 

0.2099904E— 13 

0.5609890E+05 

Ref. 48 

N 

2 

.2450268E+01 

.1066145E— 03 

— .7465337E— 07 

.1879652E— 10 

— .1025983E— 14 

.5611600E+05 

Ref. 48 

N 

3 

.2748E+01 

— .3909E— 03 

.1338E-06 

— .1191E— 10 

.3369E-15 

.5609E+05 

Ref. 47 

N 

4 

-.1227990E+01 

.1926850E— 02 

— .2437050E— 06 

.1219300E— 10 

— .1991840E— 15 

.5609000E+05 

b 

N 

5 

.1552020E+02 

— .3885790E— 02 

.3228840E— 06 

— .9605270E— 11 

.9547220E— 16 

.5609000E+05 

b 

0 

1 

.2946428E+01 

— .1638166E— 02 

.2421031E— 05 

— .1602843E— 08 

.3890696E— 12 

.2914760E+05 

Ref. 48 

0 

2 

.2542059E+01 

— .2755061E— 04 

— .3102803E— 08 

.4551067E— 11 

— .4368051E— 15 

.2923080E+05 

Ref. 48 

0 

3 

.2546E+01 

— .5952E— 04 

.2701E-07 

— .2798E— 11 

.9380E-16 

.29150E+05 

Ref. 47 

0 

4 

— .9787120E— 02 

.1244970E— 02 

— .1615440E— 06 

.8037990E— 11 

— .1262400E— 15 

.2915000E+05 

b 

o 

5 

.1642810E+02 

— .393 1 300E — 02 

.2983990E— 06 

— .8161280E— 11 

.7500430E — 16 

.2915000E+05 

b 

n 2 

1 

.3674826E+01 

— .1208150E— 02 

.2324010E— 05 

— .6321755E— 09 

— .2257725E— 12 

— .1061160E+04 

Ref. 48 

n 2 

2 

.2896319E+01 

.1515486E— 02 

— .5723527E— 06 

.9980739E— 10 

— .6522355E— 14 

— .9058620E+03 

Ref. 48 

n 2 

3 

.3727E+01 

.4684E— 03 

— .1140E— 06 

.1154E-10 

— .3293E— 15 

— .1043E+04 

Ref. 47 

n 2 

4 

.9637690E+01 

— .2572840E— 02 

.3301980E— 06 

— .1431490E— 10 

.2033260E— 15 

— .1043000E+04 

b 

n 2 

5 

— .5168080E+01 

.2333690E— 02 

— .1295340E— 06 

.2787210E— 11 

— .2135960E— 16 

— . 1043000E+04 

b 

o 2 

1 

.3625598E+01 

— .1878218E— 02 

.7055454E— 05 

— .6763513E— 08 

.2155599E— 11 

— .1047520E+04 

Ref. 48 

°2 

2 

.3621953E+01 

.7361826E— 03 

— .1965222E— 06 

.3620155E— 10 

— .2894562E— 14 

— .1201980E+04 

Ref. 48 

o 2 

3 

.3721E+01 

.4254E— 03 

-.2835E-07 

.6050E-12 

— .5186E— 17 

— .1044E+04 

Ref. 47 

o 2 

4 

.3486660E+01 

•5238420E— 03 

— .3912340E— 07 

.1009350E— 11 

— .8871830E— 17 

— .1044000E+04 

b 

o 2 

5 

.3961980E+01 

.3944550E— 03 

— .2950580E— 07 

.7397450E— 12 

— .6420930E— 17 

— .1044000E+04 

b 

NO 

1 

.4045952E+01 

— .3418178E— 02 

.7981919E— 05 

— .6113931E— 08 

.1591907E— 11 

.9745390E+04 

Ref. 48 

NO 

2 

.3189000E+01 

.1338228E— 02 

— .5289932E— 06 

.9591933E— 10 

— .6484793E— 14 

.9828330E+04 

Ref. 48 

NO 

3 

.3845E+01 

.2521E-03 

— .2658E— 07 

.2162E— 11 

— .6381E— 16 

.9764000E+04 

Ref. 47 

NO 

4 

.4330870E+01 

— .5808630E— 04 

.2805950E-07 

— .1569410E— 11 

.2410390E— 16 

.9764000E+04 

b 

NO 

5 

.2350750E+01 

.5864300E— 03 

— .3131650E— 07 

.6049510E— 12 

— .4055670E— 17 

.9764000E+04 

b 


“Ranges as follows: 1—300 < T < 1000; 2—1000 < T < 6000; 3—6000 < T < 15 000; 4—15000 < T < 25 000; 5—25000 
< T < 35 000. 

^Previously unpublished data from ref. 18. 
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Table II. Concluded 


(b) Ions and electrons 


Species 

Range 

(a) 

A i 

A 2 

a 3 

a 4 

A 5 

^6 

Source 

N+ 

1 

0.2727E+01 

— 0.2820E— 03 

0.1105E— 06 

— 0.1551E— 10 

0.7847E— 15 

0.2254E+06 

Ref. 47 

N+ 

2 

.2727E+01 

— .2820E— 03 

.1105E-06 

— .1551E— 10 

.7847E— 15 

.2254E+06 

Ref. 47 

N+ 

3 

•2499E+01 

-.3725E-05 

.1147E-07 

— .1102E— 11 

.3078E-16 

.2254E+06 

Ref. 47 

N+ 

4 

.2385610E+01 

-8349470E— 04 

— .5881510E— 08 

.1884970E— 12 

— .1611950E— 17 

.2254E+06 

6 

N+ 

5 

.2228570E+01 

.1245820E— 03 

-.8763570E-08 

.2620400E— 12 

-.2167420E-17 

.2254E+06 

b 

o+ 

1 

.2498479E+01 

.1141097E— 04 

— .2976139E— 07 

.3224653E— 10 

— .1237551E— 13 

.1879490E+06 

Ref. 48 

o+ 

2 

.2506048E+01 

— .1446424E— 04 

.1244604E— 07 

— .4685847E— 11 

.6554887E-15 

• 1879470E+06 

Ref. 48 

0+ 

3 

•2944E+01 

-.4108E-03 

.9156E-07 

— .5848E— 11 

.1190E-15 

.1879000E+06 

Ref. 47 

o+ 

4 

• 1278400E+01 

.4086590E—03 

-.2173100E-07 

.3325180E— 12 

.6316040E— 18 

.1879000E+06 

b 

o+ 

5 

.1288860E+01 

.4334250E— 03 

— .2675820E— 07 

.6215900E— 12 

-.4513150E-17 

.1879000E+06 

b 

N+ 

1 

.3397000E+01 

.4525000E— 03 

.1272000E— 06 

— .3879000E— 10 

.2459000E— 14 

.1826000E+06 

b 

N+ 

^2 

2 

•3397390E+01 

.4524870E— 03 

.1272300E— 06 

— .3879340E— 10 

.2458950E— 14 

.1826000E+06 

Ref. 47 

N+ 

3 

.3369950E+01 

.8628820E— 03 

— .1275510E— 06 

.8087120E— 11 

— .1879660E— 15 

.1826000E+06 

Ref. 47 

N+ 

^2 

4 

.4394250E+01 

.1886760E— 03 

— .7127180E— 08 

— .1751090E— 12 

.6717580E— 17 

.1826000E+06 

b 


5 

.3949290E+01 

.3679480E— 03 

— .2691020E— 07 

.6711050E— 12 

— .5824370E— 17 

.1826000E+06 

b 


1 

.3243000E+01 

.1174000E— 02 

— .3900000E— 06 

.5437000E— 10 

— .2392000E — 14 

.1400000E+06 

b 

°2 

2 

.3242980E+01 

.1173910E— 02 

— .3900420E— 06 

.5437260E— 10 

— .2392320E— 14 

.1400000E+06 

Ref. 47 


3 

.5168650E+01 

— . 86 19690E — 03 

.2041410E— 06 

— .1300410E— 10 

.2494210E— 15 

.1400000E+06 

Ref. 47 

ot 

4 

— .2801710E+ 00 

.1667410E— 02 

— .1210740E— 06 

.3211290E— 11 

— .2834890E— 16 

.1400000E+06 

b 

°2 

5 

.2044550E+01 

.1031320E— 02 

— .7404630E— 07 

.1925750E— 11 

— .1746100E— 16 

.1400000E+06 

b 

NO+ 

1 

.3668506E+01 

— .1154458E— 02 

.2175561E— 05 

-.4822747E-09 

— .2784791E— 12 

.1180340E+06 

Ref. 48 

NO+ 

2 

.2888549E+01 

.1521712E— 02 

— .5753124E— 06 

.1005108E— 09 

— .6604429E— 14 

.1181920E+06 

Ref. 48 

NO+ 

3 

.2214170E+01 

.1776060E— 02 

— .4303860E— 06 

.4173770E-10 

— .1282890E— 14 

.1181920E+06 

Ref. 48 

NO+ 

4 

— .3324050E+01 

.2441960E— 02 

— .1905720E—06 

.6858000E— 11 

— .9911240E — 16 

.1181920E+06 

b 

NO+ 

5 

-.4348760E+01 

.2401210E— 02 

— .1445990E— 06 

.3381320E— 11 

— .2825510E— 16 

.1181920E+06 

b 

e~ 

1 

.2500000E+01 

0 

0 

0 

0 

— ,7453750E-f 03 

Ref. 48 

e~ 

2 

.2500000E+01 

0 

0 

0 

0 

— .7453750E+03 

Ref. 48 

e“ 

3 

.2508E+01 

— .6332E— 05 

.1364E-08 

— .1094000E— 12 

.2934E-17 

— .7450000E+03 

Ref. 47 

e” 

4 

.250010E+01 

— .311281E— 09 

.357207E— 13 

— . 16036700E— 17 

-250707E— 22 

-.7450000E+03 

6 

e~ 

5 

.250010E+01 

.301577E— 09 

— .226204E— 13 

.667344E— 18 

— .689169E— 23 

-.7450000E+03 

b 


a Ranges as follows: 1—300 <T< 1000; 2—1000 < T < 6000; 3—6000 < T < 15000; 4—15000 <T < 25000; 5—25000 <T< 35000. 
^ Previously unpublished data from ref. 18. 
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Table III. Kinetic Model of Park 


r 

Reaction 

C f>r 

n f,r 

Ef, r/fc 


B 2 

B 3 

B 4 

B l 

1 

0 2 + M 20 + M (M = N,0) 

2.900E+23 

-2.00 

5.975E+04 

2.855 

0.988 

-6.181 

-0.023 

-0.001 

2 

0 2 + M <-► 20 + M (M = N 2 » 0 2 , NO, ions) 

9.680E+22 

-2.00 

5.975E+04 

2.855 

.988 

-6.181 

-.023 

-.001 

3 

N 2 + N <-+ 2N + N 

1.600E+22 

-1.60 

1.132E+05 

1.858 

-1.325 

-9.856 

-.174 

.008 

4 

N 2 + O 2N + O 

4.980E+22 

-1.60 

1.132E+05 

1.858 

-1.325 

-9.856 

-.174 

.008 

5 

N 2 + M «-+ 2N + M (M = N 2 , 0 2 ) 

3.700E+21 

-1.60 

1.132E+05 

1.858 

-1.325 

-9.856 

-.174 

.008 

6 

N 2 + NO «-> 2N + NO 

4.980E+21 

-1.60 

1.132E+05 

1.858 

-1.325 

-9.856 

-.174 

.008 

7 

N 2 + ions ^ 2N + ions 

8.300E+24 

-1.60 

1.132E+05 

1.858 

-1.325 

-9.856 

-.174 

.008 

8 

NO + MhN + 0 + M(M = N, O, N 2 , 0 2 , NO, ions) 

7.950E+23 

-2.00 

7.550E+04 

.792 

-.492 

-6.761 

-.091 

.004 

9 

NO + O <-> 0 2 + N 

8.370E+12 

0 

1.945E+04 

-2.063 

-1.480 

-.580 

-.114 

.005 

10 

N 2 + O «-> NO + N 

6.440E+17 

-1.00 

3.837E+04 

1.066 

-.833 

-3.095 

-.084 

.004 

11 

oj + o <-*■ o 2 + 0+ 

6.850E+13 

-.52 

1.860E+04 

-.276 

.888 

-2.180 

.055 

-.003 

12 

N 2 + N+ «-+ N+ + N 

9.850E+I2 

-.18 

1.210E+04 

.307 

-1.076 

-.878 

-.004 

-.001 

13 

NO+ + O <-► NO + 0 + 

2.750E+13 

.01 

5.100E+04 

.148 

-1.011 

-4.121 

-.132 

.006 

14 

n 2 + 0 + <-> nJ + o 

6.330E+13 

-.21 

2.220E+04 

2.979 

.382 

-3.237 

.168 

-.009 

15 

NO+ + 0 2 <-» NO + Oj 

1.030E+16 

-.17 

3.240E+04 

.424 

-1.098 

-1.941 

-.187 

.009 

16 

NO+ + N <->• N+ + O 

1.700E+13 

.40 

3.550E+04 

2.061 

.204 

-4.263 

.119 

-.006 

17 

N + O NO+ + e~ 

1.530E+09 

.37 

3.200E+04 

-7.053 

-.532 

-4.429 

.150 

-.007 

18 

O + O <-» oj + e~ 

3.850E+09 

.49 

8.060E+04 

-8.692 

-3.110 

-6.950 

-.151 

.007 

19 

N + N <-► N^" + e~ 

1.790E+09 

.77 

6.750E+04 

-4.992 

-.328 

-8.693 

.269 

-.013 

20 

O + e~ •<-> 0+ + e~ +e~ 

3.900E+33 

-3.78 

1.585E+05 

-6.113 

-2.035 

-15.311 

-.073 

.004 

21 

N + e‘«N + +«'+e- 

2.500E+33 

. 

-3.82 

1.686E+05 

-3.441 

-.577 

. 

-17.671 

.099 

-.005 
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Table IV. Kinetic Model of Dunn and Kang 


r 

Reaction 

C f,r 

n f,r 

E f,r/k 

C b , r 

n 6,r 

E b,r/k 

1 

0 2 + M <-► 20 + M (M = N, NO) 

3.600E+18 

-1.00 

5.950E+04 

3.000E+15 

-0.50 

0.000E+00 

2 

O2 + 0 20 + O 

9.000E+19 

-1.00 

5.950E+04 

7.500E+16 

-.50 

.OOOE+OO 

3 

O2 + O2 «-> 20 + O2 

3.240E+19 

-1.00 

5.950E+04 

2.700E+16 

-.50 

.000E+00 

4 

O2 + N 2 <-» 20 + N 2 

7.200E+18 

-1.00 

5.950E+04 

6.000E+15 

-.50 

•OOOE+OO 

5 

N 2 + M 2N + M (M = O, NO, 0 2 ) 

1.900E+17 

-.50 

1.130E+05 

1.100E+16 

-.50 

.000E+00 

6 

N 2 + N <-> 2N + N 

4.085E+22 

-1.50 

1.130E+05 

2.270E+21 

-1.50 

.OOOE+OO 

7 

N 2 + N 2 <-► 2N -r N 2 

4.700E+17 

-.50 

1.130E+05 

2.720E+16 

-.50 

.000E+00 

8 

NO + M<->N + 0 + M(M = 0 2 , N 2 ) 

3.900E+20 

-1.50 

7.550E+04 

1.000E+20 

-1.50 

.000E+00 

9 

NO + M«N + 0 + M(M = 0, N, NO) 

7.800E+21 

-1.50 

7.550E+04 

2.000E+21 

-1.50 

.OOOE+OO 

10 

NO + O <-*■ O2 4* N 

3.200E+09 

1.00 

1.970E+04 

1.300E+10 

1.00 

3.580E+03 

11 

N 2 + O <-> NO + N 

7.000E+13 

0 

3.800E+04 

1.560E+13 

0 

.000E+00 

12 

oj + O h o 2 + o+ 

2.920E+18 

-1.11 

2.800E+04 

7.800E+11 

.50 

.000E+00 

13 

N 2 + N+ <-> N+ + N 

2.020E+11 

.81 

1.300E+04 

7.800E+11 

.50 

•OOOE+OO 

14 

NO+ + OhNO + 0+ 

3.630E+15 

-.60 

5.080E+04 

1.500E+13 

0 

•OOOE+OO 

15 

N 2 + 0+ <- N+ + O 

3.400E+19 

-2.00 

2.300E+04 

2.480E+19 

-2.20 

•OOOE+OO 

16 

NO+ + 0 2 <-* NO + O^ 

1.800E+15 

.17 

3.300E+04 

1.800E+13 

.50 

•OOOE+OO 

17 

NO+ + N <-> NO + N+ 

1.000E+19 

-.93 

6.100E+04 

4.800E+14 

0 

•OOOE+OO 

18 

N + O «■ NO+ + e~ 

1.400E+06 

1.50 

3.190E+04 

6.700E+21 

-1.50 

•OOOE+OO 

19 

O + O <-> oj + e~ 

1.600E+17 

-.98 

8.080E+04 

8.000E+21 

-1.50 

•OOOE+OO 

20 

N + N <-> N+ + e~ 

1.400E+13 

0 

6.780E+04 

1.500E+22 

-1.50 

•OOOE+OO 

21 

O + e 0+ + e~ + e~ 

3.600E+31 

-2.91 

1.580E+05 

2.200E+40 

-4.50 

•OOOE+OO 

22 

N + e~ «-> N+ + e — + e — 

1.100E+32 

-3.14 

1.690E+05 

2.200E+40 

-4.50 

•OOOE+OO 

23 

0 2 + N 2 «-» NO + NO+ + e“ 

1.380E+20 

-1.84 

1.410E+05 

1.000E+24 

-2.50 

•OOOE+OO 

24 

N 2 + NO «N 2 + NO+ + e“ 

2.200E+15 

-.35 

1.080E+05 

2.200E+26 

-2.50 

•OOOE+OO 

25 

NO+ + O «-> 0 2 + N+ 

1.340E+13 

.31 

7.727E+04 

1.000E+14 

0 

•OOOE+OO 

26 

0 2 + NO <-> NO"*" + 0 2 + e~ 

8.800E+16 

-.35 

1.080E+05 

8.800E+26 

-2.50 

•OOOE+OO 
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Table V. Constants for Curve Fits of Electron-Neutral Energy Exchange 

Cross Section, a es 


s 

d s 

b s 

C S 

N 

5E-20 

0 

0 

0 

1.2E-20 

1.7E-24 

-2E-29 

n 2 

7.5E-20 

5.5E-24 

—IE— 28 

o 2 

2E-20 

6E-24 

0 

NO 

IE- 19 

0 

0 















Table VI. Concluded 
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z, m 


(a) Density, p/p OO ■ 

D&K T d = V TT v Case 1 

D&K T d = T Case II 



z, m 

(b) Pressure, p/pooV£>. 

Figure 1. Stagnation streamline distributions across shock layer of axisymmetric approximation to Aeroassist 
Flight Exper im ent vehicle from three different chemical kinetic models. 
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Case I 


D&K T d = ^rT v 

D&K T d = T Case II 



(a) N mole fraction. 


D&K T d = -\/T T ” Case I 

D&K T d = T Case II 



(b) O mole fraction. 

Figure 4. Profile predictions across stagnation streamline for noncontinuum, Direct-Simulation Monte Carlo 
algorithm and continuum LAURA algorithm. 
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z, m 

(a) Translational-rotational temperature. 



z, m 

(b) Vibrational-electronic temperature. 

Figure 5. Profile predictions across stagnation streamline with three different grids for case L 
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